The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
117,509 characters
Estimating Time-Varying Networks for High-Dimensional Time Series
\title{Estimating Time-Varying Networks for High-Dimensional Time Series}
\author{{\normalsize Jia Chen\thanks{Department of Economics and Related Studies,
University of York, UK. Jia Chen's research was partially supported by the
ESRC (Grant Reference: ES/T01573X/1).},\ \ \ Degui Li\thanks{Department of
Mathematics, University of York, UK. },\ \ \ Yuning Li\thanks{Department of
Economics and Related Studies, University of York, UK.},\ \ \ Oliver
Linton\thanks{Faculty of Economics, University of Cambridge, Cambridge, UK. }}\\
{\normalsize\em University of York and University of Cambridge}}
\date{{\normalsize Version: \today}}
\maketitle
\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.
\bigskip
\noindent\emph{Keywords}: CLIME, factor model, Granger causality, LASSO, local
linear smoothing, partial correlation, time-varying network, VAR.
\newpage
\section{Introduction}
\label{sec1} \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 \citep[e.g.,][]{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 \emph{static} network, which is assumed to be
invariant over time \citep[e.g.,][]{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, \cite{KSAX10} and \cite{ZLW10} study
dynamic network models with smooth time-varying structural changes; whereas
\cite{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 \citep[e.g.,][]{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
\citep[e.g.,][]{BM15, HLL15, KC15, DZZ16}. \cite{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, \cite{CFZ20}
and \cite{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,
\cite{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
\cite{SS22} simultaneously detect breaks and estimate transition matrices in
high-dimensional VAR via a three-stage procedure using the total variation
penalty. \cite{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 \cite{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. \cite{D97} is among the first to
introduce a locally stationary time series model via a time-varying spectral
representation. \cite{DS06} study a time-varying ARCH model and propose a
kernel-weighted quasi-maximum likelihood estimation method. \cite{HL10}
further consider a time-varying version of GARCH model and introduce a
semiparametric method to estimate both the parametric and nonparametric
components involved. \cite{Vo12} and \cite{ZW12} study nonparametric
kernel-based estimation and inference in a general class of locally stationary
time series. \cite{KL12} extend the locally stationary model framework to the
diffusion process. \cite{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 \citep[e.g.,][]{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 \citep[e.g.,][]{WX09, L12, FMD14, LLW14, LKZ15},
whereas the time-varying CLIME is a natural extension of the conventional
CLIME for static precision matrix estimation \citep[e.g.,][]{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
\citep[e.g.,][]{CR83} or its time-varying version \citep[e.g.,][]{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 \citep[e.g.,][]{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{sec2} introduces
the time-varying VAR and network model structures. Section \ref{sec3} presents
the estimation procedures for the time-varying transition and precision
matrices and Section \ref{sec4} gives the asymptotic properties of the
developed estimates. Section \ref{sec5} considers the factor-adjusted
time-varying VAR model and network estimation. Sections \ref{sec6} and
\ref{sec7} report simulation studies and an empirical application,
respectively. Section \ref{sec8} 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.
\section{Time-varying VAR and network models}
\label{sec2} \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{sec5} will further generalise them to the
factor-adjusted time-varying VAR and network setting.
\subsection{Time-varying VAR models}
\label{sec2.1}
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$:
\begin{equation}
X_{t}=\sum_{k=1}^{p}{\mathbf{A}}_{t,k}X_{t-k}+e_{t}\ \ \mathrm{with}
\ \ e_{t}={\boldsymbol{\Sigma}}_{t}^{1/2}\varepsilon_{t}
,\ \ t=1,\mathcal{\ldots},n,\label{eq2.1}
\end{equation}
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{eq2.1}) 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{eq2.1}) becomes the same model as that considered by
\cite{DQC17}. Furthermore, when both ${\mathbf{A}}_{t,k}$,
$k=1,\mathcal{\ldots},p$, and ${\boldsymbol{\Sigma}}_{t}$ are time-invariant
constant matrices, (\ref{eq2.1}) becomes the high-dimensional stable VAR:
\begin{equation}
X_{t}=\sum_{k=1}^{p}{\mathbf{A}}_{k}X_{t-k}+{\boldsymbol{\Sigma}}
^{1/2}\varepsilon_{t},\label{eq2.2}
\end{equation}
which has been extensively studied in the recent literature
\citep[e.g.,][]{BM15, HLL15, KC15, BB19, LZ21}. Throughout the paper, we
assume that the following conditions are satisfied.
\begin{assumption}
\label{ass:1}
\emph{(i)\ Uniformly over $\tau\in[0, 1]$, it holds that $\mathsf{det}\left(
{\mathbf{I}}_{d}-\sum_{k=1}^{p}{\mathbf{A}}_{k}(\tau)z^{k}\right) \neq0$ for
any $z\in{\mathbb{C}}$ with modulus no larger than one, where ${\mathbb{C}}$
denotes the set of complex numbers. Each entry in ${\mathbf{A}}_{k}(\cdot)$ is
second-order continuously differentiable over $[0,1]$. }
\emph{(ii)\ The precision matrix ${\boldsymbol{\Omega}}(\tau)$ is positive
definite uniformly over $\tau\in[0, 1]$, and the operator norm of
${\boldsymbol{\Sigma}}(\tau)$ is uniformly bounded over $\tau\in[0, 1]$.
Furthermore, each entry in ${\boldsymbol{\Sigma}}(\tau)$ and
${\boldsymbol{\Omega}}(\tau)$ is second-order continuously differentiable over
$[0,1]$.}
\emph{(iii)\ For any $d$-dimensional vector $u$ satisfying $\Vert u\Vert=1$,
$\mathsf{E}\left[ \exp\left\{ \iota_{1}(u^{^{\intercal}}\varepsilon_{t}
)^{2}\right\} \right] \leq C_{0}<\infty$, where $\iota_{1}$ and $C_{0}$ are
positive constants.}
\end{assumption}
The first condition in Assumption \ref{ass:1}(i) is a natural extension of the
stability assumption imposed on the constant transition matrices
\citep[e.g.,][]{Lu06}, indicating that the time-varying VAR process is locally
stationary/stable and leading to the following Wold representation
\begin{equation}
\label{eq2.3}X_{t}=\sum_{k=0}^{\infty}{\boldsymbol{\Phi}}_{t,k} e_{t-k},
\end{equation}
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,
\begin{equation}
\label{eq2.4}\max_{0\leq t\leq n}\Vert{\boldsymbol{\Phi}}_{t,k}\Vert\leq
C_{1}\rho^{k},
\end{equation}
where $C_{1}$ is a positive constant and $0<\rho<1$. A similar assumption can
be found in \cite{DQC17}. In some special model settings, (\ref{eq2.4}) may be
violated, and we refer the interested readers to the discussions in
\cite{BM15} and \cite{LZ21}. In fact, the condition (\ref{eq2.4}) 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{ass:1}(i)(ii) are common in kernel-based local estimation
method and theory. The sub-Gaussian moment condition in Assumption
\ref{ass:1}(iii) is not uncommon in the literature of high-dimensional feature
selection and covariance/precision matrix estimation \citep[e.g.,][]{W19}, and
is weaker than the Gaussian assumption frequently used in the high-dimensional
VAR literature \citep[e.g.,][]{BM15, KC15}.
\subsection{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 \cite{KSAX10} and
\cite{ZLW10} for independent network data. Following the stable network
analysis in \cite{BB19} and \cite{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 \cite{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{eq2.2}) 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{eq2.1}). 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
\begin{equation}
{\mathbb{E}}_{\tau}^{G}=\left\{ (i,j)\in{\mathbb{V}}\times{\mathbb{V}
}:\ \exists\ k\in\{1,2,\mathcal{\ldots},p\},\ a_{k,ij}(\tau)\neq0\right\}
.\label{eq2.5}
\end{equation}
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{eq2.1}). 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 \cite{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
\begin{equation}
\label{eq2.6}{\mathbb{E}}_{\tau}^{P}=\left\{ (i,j)\in{\mathbb{V}}
\times{\mathbb{V}}:\ \omega_{ij}(\tau)\neq0,\ i\neq j\right\} .
\end{equation}
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
\begin{equation}
{\mathbb{E}}^{G}=\left\{ (i,j)\in{\mathbb{V}}\times{\mathbb{V}}
:\ \exists\ k\in\{1,2,\mathcal{\ldots},p\}\ \mathrm{and}\ \tau\in
(0,1),\ a_{k,ij}(\tau)\neq0\right\} \label{eq2.7}
\end{equation}
and
\begin{equation}
{\mathbb{E}}^{P}=\left\{ (i,j)\in{\mathbb{V}}\times{\mathbb{V}}
:\ \exists\ \tau\in(0,1),\ \omega_{ij}(\tau)\neq0,\ i\neq j\right\}
.\label{eq2.8}
\end{equation}
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{sec3.4} below defines the discrete versions of the above uniform
networks and provide their estimates.
\section{Methodology}
\label{sec3} \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,
\begin{equation}
{\boldsymbol{\alpha}}_{i\bullet}(\cdot)=\left[ A_{1,i}^{^{\intercal}}
(\cdot),\mathcal{\ldots},A_{p,i}^{^{\intercal}}(\cdot)\right] ^{^{\intercal}
},\ \ {\mathbf{X}}_{t}=\left( X_{t}^{^{\intercal}},\mathcal{\ldots}
,X_{t-p+1}^{^{\intercal}}\right) ^{^{\intercal}},\label{eq3.1}
\end{equation}
and $\tau_{t}=t/n$. The time-varying VAR model (\ref{eq2.1}) can be
equivalently written as
\begin{equation}
x_{t,i}={\boldsymbol{\alpha}}_{i\bullet}^{^{\intercal}}(\tau_{t}){\mathbf{X}
}_{t-1}+e_{t,i}\ \ \mathrm{with}\ \ e_{t,i}=C_{i}^{^{\intercal}}(\tau
_{t})\varepsilon_{t},\ \ i=1,\mathcal{\ldots},d,\label{eq3.2}
\end{equation}
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
\citep[e.g.,][]{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.
\subsection{Preliminary time-varying LASSO estimation}
\label{sec3.1}
For $\tau\in(0,1)$, under the smoothness condition on the transition matrices
in Assumption \ref{ass:1}(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
\citep[e.g.,][]{FG96}:
\begin{equation}
\mathcal{L}_{i}({\boldsymbol{\alpha}},{\boldsymbol{\beta}}\ |\ \tau)=\frac
{1}{n}\sum\limits_{t=1}^{n}\left\{ x_{t,i}-\left[ {\boldsymbol{\alpha}
}+{\boldsymbol{\beta}}(\tau_{t}-\tau)\right] ^{^{\intercal}}{\mathbf{X}
}_{t-1}\right\} ^{2}K_{h}(\tau_{t}-\tau),\label{eq3.3}
\end{equation}
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$ \citep[e.g.,][]{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{eq3.3}).
The LASSO estimation was first introduced by \cite{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
\begin{equation}
\label{eq3.4}\mathcal{L}_{i}^{\ast}({\boldsymbol{\alpha}}, {\boldsymbol{\beta
}}\ |\ \tau)=\mathcal{L}_{i}({\boldsymbol{\alpha}}, {\boldsymbol{\beta}
}\ |\ \tau)+\lambda_{1} \left( \vert{\boldsymbol{\alpha}}\vert_{1}
+h\vert{\boldsymbol{\beta}}\vert_{1}\right) ,
\end{equation}
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.
\subsection{Penalised local linear estimation with weighted group LASSO}
\label{sec3.2}
In order to estimate the uniform Granger causality network, we next introduce
a \emph{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:
\begin{equation}
\mathcal{Q}_{i}({\mathbf{A}},{\mathbf{B}})=\sum_{t=1}^{n}\mathcal{L}
_{i}({\boldsymbol{\alpha}}_{\bullet t},{\boldsymbol{\beta}}_{\bullet
t}\ |\ \tau_{t})+\sum_{j=1}^{pd}p_{\lambda_{2}}^{\prime}\left( \left\Vert
\widetilde{\boldsymbol{\alpha}}_{i,j}\right\Vert \right) \Vert
{\boldsymbol{\alpha}}_{j}\Vert+\sum_{j=1}^{pd}p_{\lambda_{2}}^{\prime}\left(
\widetilde{D}_{i,j}\right) \Vert h{\boldsymbol{\beta}}_{j}\Vert,\label{eq3.5}
\end{equation}
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 \cite{FL01} and $I(\cdot)$ being the
indicator function. The penalty terms in (\ref{eq3.5}) are motivated by the
local linear approximation to the SCAD penalty function \citep{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{eq3.5}) serve as the
weights for the group LASSO, and their values are determined by the
preliminary estimates in Section \ref{sec3.1}, 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
\begin{align}
& \widehat{\mathbf{A}}_{i}=\left( \widehat{\boldsymbol{\alpha}}
_{i,1},\mathcal{\ldots},\widehat{\boldsymbol{\alpha}}_{i,pd}\right)
\ \ \mathrm{with}\ \ \widehat{\boldsymbol{\alpha}}_{i,j}=\left[
\widehat{\alpha}_{i,j}(\tau_{1}),\mathcal{\ldots},\widehat{\alpha}_{i,j}
(\tau_{n})\right] ^{^{\intercal}},\nonumber\\
& \widehat{\mathbf{B}}_{i}=\left( \widehat{\boldsymbol{\alpha}}_{i,1}
^{\prime},\mathcal{\ldots},\widehat{\boldsymbol{\alpha}}_{i,pd}^{\prime
}\right) \ \ \mathrm{with}\ \ \widehat{\boldsymbol{\alpha}}_{i,j}^{\prime
}=\left[ \widehat{\alpha}_{i,j}^{\prime}(\tau_{1}),\mathcal{\ldots}
,\widehat{\alpha}_{i,j}^{\prime}(\tau_{n})\right] ^{^{\intercal}}.\nonumber
\end{align}
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 \cite{LKZ15} and \cite{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.
\subsection{Estimation of the time-varying precision matrix}
\label{sec3.3}
In this section, we study the estimation of ${\boldsymbol{\Omega}}(\cdot)$ in
model (\ref{eq2.1}), 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
\citep[e.g.,][]{LF09, Y10, CLL11}. \cite{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{sec3.2}, 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
\begin{equation}
\widehat{e}_{t}=\left( \widehat{e}_{t,1},\mathcal{\ldots},\widehat{e}
_{t,d}\right) ^{^{\intercal}}=X_{t}-\sum_{k=1}^{p}\widehat{\mathbf{A}}
_{k}(\tau_{t})X_{t-k},\ \ t=1,\mathcal{\ldots},n.\label{eq3.6}
\end{equation}
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 {\small
\begin{equation}
\mathcal{S}(q,\xi_{d})=\left\{ {\mathbf{W}}(\tau)=\left[ w_{ij}
(\tau)\right] _{d\times d},0\leq\tau\leq1:\ {\mathbf{W}}(\tau)\succ
0,\ \sup_{0\leq\tau\leq1}\Vert{\mathbf{W}}(\tau)\Vert_{1}\leq C_{2}
,\ \sup_{0\leq\tau\leq1}\max_{1\leq i\leq d}\sum_{j=1}^{d}|w_{ij}(\tau
)|^{q}\leq\xi_{d}\right\} ,\label{eq3.7}
\end{equation}
} 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
\begin{equation}
\widehat{\boldsymbol{\Sigma}}(\tau)=\left[ \widehat{\sigma}_{ij}
(\tau)\right] _{d\times d}\ \ \mathrm{with}\ \ \widehat{\sigma}_{ij}
(\tau)=\sum_{t=1}^{n}\varpi_{n,t}(\tau)\widehat{e}_{t,i}\widehat{e}_{t,j}
/\sum_{t=1}^{n}\varpi_{n,t}(\tau),\label{eq3.8}
\end{equation}
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{eq3.7}), we estimate ${\boldsymbol{\Omega}}(\tau)$
via the time-varying CLIME method:
\begin{equation}
\widetilde{\boldsymbol{\Omega}}(\tau)=\left[ \widetilde{\omega}_{ij}
(\tau)\right] _{d\times d}=\operatorname*{arg\,min}_{\boldsymbol{\Omega}}|{\boldsymbol{\Omega
}}|_{1}\ \ \ \ \ \mathrm{subject\ to}\ \ \left\Vert \widehat
{\boldsymbol{\Sigma}}(\tau){\boldsymbol{\Omega}}-{\mathbf{I}}_{d}\right\Vert
_{\mathrm{max}}\leq\lambda_{3},\label{eq3.9}
\end{equation}
where $\lambda_{3}$ is a tuning parameter. As the underlying time-varying
precision matrix is symmetric, the matrix estimate obtained from (\ref{eq3.9})
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
\begin{equation}
\widehat{\omega}_{ij}(\tau)=\widehat{\omega}_{ji}(\tau)=\widetilde{\omega
}_{ij}(\tau)I\left( |\widetilde{\omega}_{ij}(\tau)|\leq|\widetilde{\omega
}_{ji}(\tau)|\right) +\widetilde{\omega}_{ji}(\tau)I\left( |\widetilde
{\omega}_{ij}(\tau)|>|\widetilde{\omega}_{ji}(\tau)|\right) .\label{eq3.10}
\end{equation}
\subsection{Estimation of uniform time-varying networks}
\label{sec3.4}
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:
\begin{equation}
{\mathbb{E}}_{n}^{G}=\left\{ (i,j)\in{\mathbb{V}}\times{\mathbb{V}}
:\ \exists\ k\in\{1,2,\mathcal{\ldots},p\}\ \ \mathrm{and}\ \ t\in
\{1,\mathcal{\ldots},n\},\ a_{k,ij}(\tau_{t})\neq0\right\} \label{eq3.11}
\end{equation}
and
\begin{equation}
{\mathbb{E}}_{n}^{P}=\left\{ (i,j)\in{\mathbb{V}}\times{\mathbb{V}}
:\ \exists\ t\in\{1,\mathcal{\ldots},n\},\ \omega_{ij}(\tau_{t})\neq0,\ i\neq
j\right\} .\label{eq3.12}
\end{equation}
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{sec3.2} and
\ref{sec3.3}, we can estimate ${\mathbb{E}}_{n}^{G}$ by
\begin{equation}
\widehat{\mathbb{E}}_{n}^{G}=\left\{ (i,j)\in{\mathbb{V}}\times{\mathbb{V}
}:\ \exists\ k\in\{1,2,\mathcal{\ldots},p\},\ \sum_{t=1}^{n}\widehat{a}
_{k,ij}^{2}(\tau_{t})>0\right\} ,\label{eq3.13}
\end{equation}
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
\begin{equation}
\widehat{\mathbb{E}}_{n}^{P}=\left\{ (i,j)\in{\mathbb{V}}\times{\mathbb{V}
}:\ \exists\ t\in\{1,\mathcal{\ldots},n\},\ \left\vert \widehat{\omega}
_{ij}(\tau_{t})\right\vert \geq\lambda_{3},\ i\neq j\right\} ,\label{eq3.14}
\end{equation}
where $\lambda_{3}$ is the tuning parameter used in the time-varying CLIME.
\section{Main theoretical results}
\label{sec4} \setcounter{equation}{0}
To ease the notational burden, throughout this section, we focus on the
time-varying VAR(1) model:
\begin{equation}
X_{t}={\mathbf{A}}(\tau_{t})X_{t-1}+{\boldsymbol{\Sigma}}_{t}^{1/2}
\varepsilon_{t},\label{eq4.1}
\end{equation}
where ${\mathbf{A}}(\tau)=\left[ \alpha_{ij}(\tau)\right] _{d\times d}$. For
a general time-varying VAR($p$) model (\ref{eq2.1}), 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{eq3.1}), ${\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(
\begin{array}
[c]{ccccc}
{\mathbf{A}}_{t,1} & {\mathbf{A}}_{t,2} & \mathcal{\ldots} & {\mathbf{A}
}_{t,p-1} & {\mathbf{A}}_{t,p}\\
{\mathbf{I}}_{d} & {\mathbf{O}}_{d\times d} & \mathcal{\ldots} & {\mathbf{O}
}_{d\times d} & {\mathbf{O}}_{d\times d}\\
\vdots & \vdots & \vdots & \vdots & \vdots\\
{\mathbf{O}}_{d\times d} & {\mathbf{O}}_{d\times d} & \mathcal{\ldots} &
{\mathbf{I}}_{d} & {\mathbf{O}}_{d\times d}
\end{array}
\right) .
\]
\subsection{Uniform consistency of the time-varying LASSO estimates}
\label{sec4.1}
Define {\small
\begin{equation}
\label{eq4.2}{\boldsymbol{\Psi}}(\tau)=\left[
\begin{array}
[c]{cc}
{\boldsymbol{\Psi}}_{0}(\tau) & {\boldsymbol{\Psi}}_{1}(\tau)\\
{\boldsymbol{\Psi}}_{1}(\tau) & {\boldsymbol{\Psi}}_{2}(\tau)
\end{array}
\right] \ \ \mathrm{with} \ \ {\boldsymbol{\Psi}}_{k}(\tau)=\frac{1}{n}
\sum\limits_{t=1}^{n}\left( \frac{\tau_{t}-\tau}{h}\right) ^{k} X_{t-1}
X_{t-1}^{^{\intercal}}K_{h}(\tau_{t} - \tau),\ \ k=0,1,2,
\end{equation}
} 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{sec3.2} but with $p=1$. To derive the uniform consistency
property of the preliminary time-varying LASSO estimates defined in Section
\ref{sec3.1}, we need the following assumptions, some of which may be weakened
at the cost of lengthier proofs.
\begin{assumption}
\label{ass:2}
\emph{(i)\ The kernel $K(\cdot)$ is a bounded, continuous and symmetric
probability density function with a compact support $[-1,1]$.}
\emph{(ii)\ The bandwidth $h$ satisfies
\[
nh/\log^{2} (n\vee d)\rightarrow\infty\ \ \mbox{and}\ \ sh^{2}\log(n\vee
d)\rightarrow0,
\]
where $s=\max_{1\leq i\leq d}s_{i}$ with $s_{i}$ being the cardinality of the
index set $\mathcal{J}_{i}$. }
\end{assumption}
\begin{assumption}
\label{ass:3}
\emph{(i)\ The tuning parameter $\lambda_{1}$ satisfies }
\[
\zeta_{n,d}:=\log(n\vee d)\left[ (nh)^{-1/2}+sh^{2}\right] =o(\lambda
_{1})\ \ \mbox{and}\ \ \sqrt{s}\lambda_{1}/h\rightarrow0.
\]
\emph{(ii)\ There exists a positive constant $\kappa_{0}$ such that, with
probability approaching one (w.p.a.1),}
\begin{equation}
\label{eq4.3}\min_{1\leq i\leq d}\min_{1\leq t\leq n}\inf_{u\in\mathcal{B}
_{i}(\tau_{t})}u^{^{\intercal}} {\boldsymbol{\Psi}}(\tau_{t})u\geq\kappa_{0}.
\end{equation}
\end{assumption}
Assumption \ref{ass:2}(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{ass:2}(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{ass:3}(ii) can be seen as a uniform version of the so-called restricted
eigenvalue condition widely used in high-dimensional linear regression models
\citep[e.g.,][]{BRT09,BM15}. Appendix D in the supplement provides sufficient
conditions for the high-dimensional locally stationary Gaussian time series to
satisfy Assumption \ref{ass:3}(ii). Furthermore, with the Hanson-Wright
inequality for time-varying (non-Gaussian) VAR processes
\citep[e.g., Proposition 6.2 in][]{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{eq4.3}) 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}
\begin{theorem}\label{thm:4.1}
Suppose that Assumptions \ref{ass:1}--\ref{ass:3} are satisfied. Then we have
\begin{equation}\label{eq4.4}
\max_{1\leq i\leq d}\max_{1\leq t\leq n}\left\Vert \widetilde{\boldsymbol\alpha}_{i\bullet}(\tau_t)-{\boldsymbol\alpha}_{i\bullet}(\tau_t)\right\Vert =O_P\left(\sqrt{s}\lambda_1\right).
\end{equation}
\end{theorem}
Theorem \ref{thm:4.1} 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{eq4.4}) can be strengthened to
\begin{equation}
\label{eq4.5}\max_{1\leq i\leq d}\sup_{0\leq\tau\leq1}\left\Vert
\widetilde{\boldsymbol{\alpha}}_{i\bullet}(\tau)-{\boldsymbol{\alpha}
}_{i\bullet}(\tau)\right\Vert =O_{P}\left( \sqrt{s}\lambda_{1}\right) .
\end{equation}
A similar uniform convergence property holds for the first-order derivative
function estimates, see (A.1) in the proof of Theorem \ref{thm:4.1}.
\subsection{The oracle property of the weighted group LASSO estimates}
\label{sec4.2}
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:
\begin{align}
& \widehat{\mathbf{A}}_{i}^{o}=\left( \widehat{\boldsymbol{\alpha}}_{i,1}
^{o},\mathcal{\ldots},\widehat{\boldsymbol{\alpha}}_{i,d}^{o}\right)
\ \ \mathrm{with}\ \ \widehat{\boldsymbol{\alpha}}_{i,j}^{o}=\left[
\widehat{\alpha}_{i,j}^{o}(\tau_{1}),\mathcal{\ldots},\widehat{\alpha}
_{i,j}^{o}(\tau_{n})\right] ^{^{\intercal}},\label{eq4.6}\\
& \widehat{\mathbf{B}}_{i}^{o}=\left( \widehat{\boldsymbol{\alpha}}
_{i,1}^{\prime o},\mathcal{\ldots},\widehat{\boldsymbol{\alpha}}_{i,d}^{\prime
o}\right) \ \ \mathrm{with}\ \ \widehat{\boldsymbol{\alpha}}_{i,j}^{\prime
o}=\left[ \widehat{\alpha}_{i,j}^{\prime o}(\tau_{1}),\mathcal{\ldots
},\widehat{\alpha}_{i,j}^{\prime o}(\tau_{n})\right] ^{^{\intercal}
},\label{eq4.7}
\end{align}
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.
\begin{assumption}
\label{ass:4}
\emph{(i)\ The tuning parameter $\lambda_{2}$ satisfies
\[
\sqrt{n}s\log(n\vee d)\zeta_{n,d}+\sqrt{ns}\lambda_{1}=o(\lambda_{2}),
\]
where $\zeta_{n,d}$ is defined in Assumption \ref{ass:3}(i).}
\emph{(ii)\ It holds that
\[
\min_{1\leq i\leq d}\min_{j\in\mathcal{J}_{i}}\left( \sum_{t=1}^{n}
\alpha_{i,j}^{2}(\tau_{t})\right) ^{\frac{1}{2}}\geq(a_{0}+1)\lambda
_{2}\ \ \mbox{and}\ \ \min_{1\leq i\leq d}\min_{j\in\mathcal{J}_{i}^{\prime}
}D_{i,j}\geq(a_{0}+1)\lambda_{2},
\]
where $a_{0}=3.7$ is defined in the SCAD penalty.}
\end{assumption}
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{ass:4}(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{ass:4}(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 \cite{LKZ15} and Assumption 8 in \cite{CLWZ21}, and they
ensure that the significant coefficient functions and derivatives can be
detected \emph{w.p.a.1}.
\begin{theorem}\label{thm:4.2}
Suppose that Assumptions \ref{ass:1}--\ref{ass:4} are satisfied. The minimiser to the objective function of the weighted group LASSO, ${\cal Q}_{i}({\mathbf A}, {\mathbf B})$, exists and equals the oracle estimates defined in (\ref{eq4.6}) and (\ref{eq4.7}) w.p.a.1. In addition, we have the following mean squared convergence result:
\begin{equation}\label{eq4.8}
\max_{1\leq i\leq d}\frac{1}{n}\sum_{t=1}^n\sum_{j=1}^d\left[ \widehat{\alpha}_{ij}(\tau_t)-\alpha_{ij}(\tau_t)\right]^2=O_P\left(s\zeta_{n,d}^2\right),
\end{equation}
where $s$ is defined in Assumption \ref{ass:2}(ii) and $\zeta_{n,d}$ is defined in Assumption \ref{ass:3}(i).
\end{theorem}
\medskip
Since the penalised local linear estimates are identical to the infeasible
oracle estimates defined in (\ref{eq4.6}) and (\ref{eq4.7}) \emph{w.p.a.1},
the sparsity property holds for the global model selection procedures proposed
in Section \ref{sec3.2}, i.e., the zero elements in the time-varying
transition matrix can be estimated exactly as zeros. Following the proof of
Theorem \ref{thm:4.2}, we may verify properties (i)--(iv) for the folded
concave penalty function discussed in \cite{FXZ14} \emph{w.p.a.1}. Hence,
Theorem \ref{thm:4.2} may be regarded as a generalisation of Theorem 1 in
\cite{FXZ14} and Theorem 3.1 in \cite{LKZ15} to high-dimensional locally
stationary time series.
\medskip
With the oracle property in Theorem \ref{thm:4.2}, it is straightforward to
derive the following consistency property of the network estimates for the
directed edges of Granger causality linkages.
\setcounter{corollary}{0}
\begin{corollary}
\label{cor:4.1}
\emph{Under the assumptions of Theorem \ref{thm:4.2}, we have}
\begin{equation}
\label{eq4.9}\mathsf{P}\left( \widehat{\mathbb{E}}_{n}^{G}={\mathbb{E}}
_{n}^{G}\right) \rightarrow1.
\end{equation}
\end{corollary}
\subsection{Uniform consistency of the time-varying CLIME estimates}
\label{sec4.3}
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}$.
\begin{assumption}
\label{ass:5}
\emph{(i)\ The bandwidth $b$ satisfies
\[
b\rightarrow0\ \ \ \mbox{and}\ \ \ nb/[\log(n\vee d)]^{3}\rightarrow\infty.
\]
In addition, $s\zeta_{n,d}\sqrt{\log(n\vee d)}\rightarrow0$, where
$\zeta_{n,d}$ is defined in Assumption \ref{ass:3}(i).}
\emph{(ii)\ There exists a sufficiently large constant $C_{3}$ such that
$\lambda_{3}=C_{3}\left( \nu_{n,d}^{\diamond}+\nu_{n,d}^{\ast}\right) $, where
}
\[
\nu_{n,d}^{\diamond}=\left[ \frac{\log(n\vee d)}{nb}\right] ^{1/2}
+b^{2}\ \ \ \mbox{and}\ \ \ \nu_{n,d}^{\ast}=s\zeta_{n,d}\sqrt{\log(n\vee
d)}.
\]
\end{assumption}
The following theorem gives the uniform convergence rates of the time-varying
precision matrix estimate $\widehat{\boldsymbol{\Omega}}(\tau)$ under various
matrix norms.
\begin{theorem}\label{thm:4.3}
Suppose Assumptions \ref{ass:1}--\ref{ass:5} are satisfied and $\left\{{\boldsymbol\Omega}(\tau): 0\leq \tau\leq 1\right\}\in{\cal S}(q, \xi_d)$. Then we have
\begin{eqnarray}
&&\sup_{0\leq \tau\leq 1}\left\Vert\widehat{\boldsymbol\Omega}(\tau)-{\boldsymbol\Omega}(\tau)\right\Vert_{\max}=O_P\left(\nu_{n,d}^\diamond+\nu_{n,d}^\ast\right),\label{eq4.10}\\
&&\sup_{0\leq \tau\leq 1}\left\Vert \widehat{\boldsymbol\Omega}(\tau)-{\boldsymbol\Omega}(\tau)\right\Vert =O_P\left( \xi_d(\nu_{n,d}^\diamond+\nu_{n,d}^\ast)^{1-q}\right),\label{eq4.11}\\
&&\sup_{0\leq\tau\leq1}\frac{1}{d}\left\Vert \widehat{\boldsymbol\Omega}(\tau)-{\boldsymbol\Omega}(\tau)\right\Vert _{F}^2=O_P\left( \xi_d(\nu_{n,d}^\diamond+\nu_{n,d}^\ast)^{2-q}\right),\label{eq4.12}
\end{eqnarray}
where $\xi_d$ is defined in (\ref{eq3.7}), $\nu_{n,d}^\diamond$ and $\nu_{n,d}^\ast$ are defined in Assumption \ref{ass:5}(ii).
\end{theorem}
\medskip
The uniform convergence rates in Theorem \ref{thm:4.3} 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{thm:4.3} can be seen as an extension of Theorem 1 in
\cite{CLL11} to the high-dimensional locally stationary time series setting.
\medskip
From Theorem \ref{thm:4.3}, we readily have the following consistency property
for the network estimates of the undirected edges of partial correlation linkages.
\begin{corollary}
\label{cor:4.2}
\emph{Under the assumptions of Theorem \ref{thm:4.3}, if $\min_{(i,j)\in
{\mathbb{E}}^{P}} \min_{1\leq t\leq n}\vert\omega_{ij}(\tau_{t})\vert
\gg\lambda_{3}$, we have}
\begin{equation}
\label{eq4.13}\mathsf{P}\left( \widehat{\mathbb{E}}_{n}^{P}={\mathbb{E}}
_{n}^{P}\right) \rightarrow1.
\end{equation}
\end{corollary}
\section{Factor-adjusted time-varying VAR and networks}
\label{sec5} \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:
\begin{equation}
Z_{t}={\boldsymbol{\Lambda}}F_{t}+X_{t},\ \ t=1,\mathcal{\ldots}
,n,\label{eq5.1}
\end{equation}
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{eq2.1}). More generally, we may assume the
following time-varying factor model structure:
\begin{equation}
Z_{t}={\boldsymbol{\Lambda}}_{t}F_{t}+X_{t},\ \ t=1,\mathcal{\ldots
},n,\label{eq5.2}
\end{equation}
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
\citep[e.g.,][]{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{sec3}. 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 \citep[e.g.,][]{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{eq5.1}), 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
\begin{equation}
\widehat{X}_{t}=Z_{t}-\widehat{\boldsymbol{\Lambda}}\widehat{F}_{t}
,\ \ t=1,\mathcal{\ldots},n.\label{eq5.3}
\end{equation}
For the time-varying factor model (\ref{eq5.2}), 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{sec3.1}, 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
\begin{equation}
\widehat{X}_{t}=Z_{t}-\widehat{\boldsymbol{\Lambda}}(\tau_{t})\widehat{F}
(\tau_{t}),\ \ t=1,\mathcal{\ldots},n,\label{eq5.4}
\end{equation}
where we've kept the same notation $\widehat{X}_{t}$ as in (\ref{eq5.3}) to
avoid notational burden.
As in Section \ref{sec4}, 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{sec3}.
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{eq3.13}) and (\ref{eq3.14}), 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{ass:3}--\ref{ass:5} to incorporate the approximation error of the
idiosyncratic error components.
\begin{assumption}
\label{ass:6}
\emph{(i) Denote $\delta_{X}=\max_{1\leq t\leq n}\left\vert \widehat{X}
_{t}-X_{t}\right\vert _{\max}$. It holds that $[\log(n\vee d)]^{1/2}
s\delta_{X}=o_{P}(1)$.}
\emph{(ii) Assumption \ref{ass:3}(i) holds when $\zeta_{n,d}$ is replaced by
$\zeta_{n,d}^{\dagger}=\zeta_{n,d}+[\log(n\vee d)]^{1/2}s\delta_{X}$.}
\emph{(iii) Assumption \ref{ass:4}(i) holds when $\zeta_{n,d}$ is replaced by
$\zeta_{n,d}^{\dagger}$.}
\emph{(iv) Assumption \ref{ass:5} holds when $\zeta_{n,d}$ and $\nu
_{n,d}^{\ast}$ are replaced by $\zeta_{n,d}^{\dagger}$ and $\nu_{n,d}
^{\dagger}=s\zeta_{n,d}^{\dagger}\sqrt{\log(n\vee d)}$, respectively.}
\end{assumption}
Assumption \ref{ass:6}(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 \cite{FLM13}, a typical rate for
the approximation error from PCA estimation of the conventional factor model
(\ref{eq5.1}) is
\begin{equation}
\label{eq5.5}\delta_{X}=O_{P}\left( (\log n)^{1/2}\left[ (\log d)^{1/2}
n^{-1/2}+n^{1/\upsilon}d^{-1/2}\right] \right) ,
\end{equation}
where $\upsilon>2$ is a positive number related to moment restrictions. From
Theorem 3.5 in \cite{SW17}, we may obtain the typical uniform rate for
$\delta_{X}$ under the time-varying factor model (\ref{eq5.2}) when the local
PCA estimation is used. In Assumption \ref{ass:6}(ii)--(iv), we amend
Assumptions \ref{ass:3}(i), \ref{ass:4}(i) and \ref{ass:5}(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{eq5.5}) 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{ass:6}(ii)--(iv) and $\delta_{X}$ would not be involved in the
estimation convergence rates under model (\ref{eq5.1}).
The following two propositions extend the theoretical results in Section
\ref{sec4} to the factor-adjusted time-varying VAR and networks.
\setcounter{prop}{0}
\begin{prop}
\label{prop:5.1}
\emph{Suppose that the factor model (\ref{eq5.1}) or (\ref{eq5.2}), and
Assumptions \ref{ass:1}, \ref{ass:2} and \ref{ass:3}(ii) are satisfied.}
\emph{(i) Under Assumption \ref{ass:6}(i)(ii), we have}
\begin{equation}
\label{eq5.6}\max_{1\leq i\leq d}\max_{1\leq t\leq n}\sum_{j=1}^{d}\left[
\widetilde{\alpha}_{ij}^{\dagger}(\tau_{t})-{\alpha}_{ij}(\tau_{t})\right]
^{2}=O_{P}\left( s\lambda_{1}^{2}\right) .
\end{equation}
\emph{(ii) Under Assumption \ref{ass:6}(i)--(iii), the oracle property holds
for the second-stage weighted group LASSO estimates and furthermore, }
\begin{equation}
\label{eq5.7}\max_{1\leq i\leq d}\frac{1}{n}\sum_{t=1}^{n}\sum_{j=1}
^{d}\left[ \widehat{\alpha}_{ij}^{\dagger}(\tau_{t})-\alpha_{ij}(\tau
_{t})\right] ^{2}=O_{P}\left( s\left( \zeta_{n,d}^{\dagger}\right) ^{2}\right)
.
\end{equation}
\emph{(iii) Under Assumption \ref{ass:6} and the sparsity condition that
$\left\{ {\boldsymbol{\Omega}}(\tau): 0\leq\tau\leq1\right\} \in\mathcal{S}(q,
\xi_{d})$, we have}
\begin{align}
& \sup_{0\leq\tau\leq1}\left\Vert \widehat{\boldsymbol{\Omega}}^{\dagger}
(\tau)-{\boldsymbol{\Omega}}(\tau)\right\Vert _{\max}=O_{P}\left( \nu
_{n,d}^{\diamond}+\nu_{n,d}^{\dagger}\right) ,\label{eq5.8}\\
& \sup_{0\leq\tau\leq1}\left\Vert \widehat{\boldsymbol{\Omega}}^{\dagger}
(\tau)-{\boldsymbol{\Omega}}(\tau)\right\Vert =O_{P}\left( \xi_{d}(\nu
_{n,d}^{\diamond}+\nu_{n,d}^{\dagger})^{1-q}\right) ,\label{eq5.9}\\
& \sup_{0\leq\tau\leq1}\frac{1}{d}\left\Vert \widehat{\boldsymbol{\Omega}
}^{\dagger}(\tau)-{\boldsymbol{\Omega}}(\tau)\right\Vert _{F}^{2}=O_{P}\left(
\xi_{d}(\nu_{n,d}^{\diamond}+\nu_{n,d}^{\dagger})^{2-q}\right) .\label{eq5.10}
\end{align}
\end{prop}
\begin{prop}
\label{prop:5.2}
\emph{(i) Under the assumptions of Proposition \ref{prop:5.1}(ii), we have}
\begin{equation}
\label{eq5.11}\mathsf{P}\left( \widehat{\mathbb{E}}_{n}^{G,\dagger
}={\mathbb{E}}_{n}^{G}\right) \rightarrow1.
\end{equation}
\emph{(ii) Under the assumptions of Proposition \ref{prop:5.1}(iii) and
$\min_{(i,j)\in{\mathbb{E}}^{P}} \min_{1\leq t\leq n}\vert\omega_{ij}(\tau
_{t})\vert\gg\lambda_{3}$, we have}
\begin{equation}
\label{eq5.12}\mathsf{P}\left( \widehat{\mathbb{E}}_{n}^{P,\dagger
}={\mathbb{E}}_{n}^{P}\right) \rightarrow1.
\end{equation}
\end{prop}
\section{Monte-Carlo simulation}
\label{sec6} \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 ``\textsf{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 \cite{LKZ15}. The bandwidth
for the local PCA is set as $h_{\ast}=(2.35/\sqrt{12})[\sqrt{d}/n] ^{1/5}$ as
in \cite{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{eq3.3}) 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.
\medskip
\noindent\textbf{Example 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{tab:1} 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{tab:2} 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.
\begin{table}
\caption{Transition matrix and Granger network estimation in Example 1.}
\label{tab:1}
\centering
\begin{tabular}
[c]{llllllllll}\hline\hline
& & \multicolumn{2}{l}{tv-wgLASSO} & & \multicolumn{2}{l}{tv-Oracle} & &
\multicolumn{2}{l}{tv-Full}\\\cline{3-4}\cline{6-7}\cline{9-10}
measure & dimension & $n=200$ & $n=400$ & & $n=200$ & $n=400$ & & $n=200$ &
$n=400$\\\hline
FP & $d=50$ & 0.97 & 0.04 & & 0 & 0 & & 2450 & 2450\\
& $d=100$ & 1.73 & 0.08 & & 0 & 0 & & - & 9900\\
FN & $d=50$ & 3.53 & 0.08 & & 0 & 0 & & 0 & 0\\
& $d=100$ & 8.55 & 0.15 & & 0 & 0 & & - & 0\\
TPR & $d=50$ & 0.929 & 0.998 & & 1 & 1 & & 1 & 1\\
& $d=100$ & 0.915 & 0.999 & & 1 & 1 & & - & 1\\
TNR & $d=50$ & 1.000 & 1.000 & & 1 & 1 & & 0 & 0\\
& $d=100$ & 1.000 & 1.000 & & 1 & 1 & & - & 0\\
PPV & $d=50$ & 0.980 & 0.999 & & 1 & 1 & & 0.02 & 0.02\\
& $d=100$ & 0.982 & 0.999 & & 1 & 1 & & - & 0.01\\
NPV & $d=50$ & 0.999 & 1.000 & & 1 & 1 & & 1 & 1\\
& $d=100$ & 0.999 & 1.000 & & 1 & 1 & & - & 1\\
F1 & $d=50$ & 0.953 & 0.999 & & 1 & 1 & & 0.039 & 0.039\\
& $d=100$ & 0.947 & 0.999 & & 1 & 1 & & - & 0.020\\
MCC & $d=50$ & 0.953 & 0.999 & & 1 & 1 & & 0 & 0\\
& $d=100$ & 0.947 & 0.999 & & 1 & 1 & & - & 0\\\hline
average $R^{2}$ & $d=50$ & 0.289 & 0.296 & & 0.296 & 0.297 & & 0.933 &
0.721\\
& $d=100$ & 0.296 & 0.306 & & 0.305 & 0.307 & & - & 0.959\\\hline
EE$_{A}$ & $d=50$ & 0.214 & 0.160 & & 0.185 & 0.163 & & 54.29 & 1.410\\
& $d=100$ & 0.224 & 0.163 & & 0.189 & 0.166 & & - & 112.8\\
RMSE$_{e}$ & $d=50$ & 0.203 & 0.115 & & 0.162 & 0.120 & & 1.119 & 0.876\\
& $d=100$ & 0.213 & 0.113 & & 0.159 & 0.119 & & - & 1.145\\\hline
\end{tabular}
\par
\begin{flushleft}
\emph{In all the tables, except for exact values of 0's and 1's, the FP and FN
measures are rounded to 2 decimal places, while the others are rounded to 3
decimal places. }
\end{flushleft}
\end{table}
\begin{table}
\caption{Precision matrix and partial correlation network estimation in
Example 1.}
\label{tab:2}
\centering
\begin{tabular}
[c]{llllllllll}\hline\hline
& & \multicolumn{2}{c}{tv-CLIME} & & \multicolumn{2}{c}{infeasible tv-CLIME}
& & \multicolumn{2}{c}{tv-GLASSO}\\\cline{3-4}\cline{6-7}\cline{9-10}
measure & dimension & $n=200$ & $n=400$ & & $n=200$ & $n=400$ & & $n=200$ &
$n=400$\\\hline
FP & $d=50$ & 0 & 0.02 & & 0 & 0.02 & & 0 & 0\\
& $d=100$ & 0 & 0.03 & & 0 & 0.01 & & 0 & 0\\
FN & $d=50$ & 5.06 & 0 & & 3.49 & 0 & & 9.24 & 0\\
& $d=100$ & 13.25 & 0 & & 9.01 & 0 & & 28.31 & 0\\
TPR & $d=50$ & 0.798 & 1 & & 0.860 & 1 & & 0.630 & 0\\
& $d=100$ & 0.735 & 1 & & 0.820 & 1 & & 0.434 & 0\\
TNR & $d=50$ & 1 & 1.000 & & 1 & 1.000 & & 1 & 1\\
& $d=100$ & 1 & 1.000 & & 1 & 1.000 & & 1 & 1\\
PPV & $d=50$ & 1 & 0.999 & & 1 & 0.999 & & 1 & 1\\
& $d=100$ & 1 & 0.999 & & 1 & 1.000 & & 1 & 1\\
NPV & $d=50$ & 0.996 & 1 & & 0.097 & 1 & & 0.992 & 1\\
& $d=100$ & 0.997 & 1 & & 0.998 & 1 & & 0.994 & 1\\
F1 & $d=50$ & 0.884 & 1.000 & & 0.922 & 1.000 & & 0.768 & 1\\
& $d=100$ & 0.845 & 1.000 & & 0.899 & 1.000 & & 0.600 & 1\\
MCC & $d=50$ & 0.889 & 1.000 & & 0.925 & 1.000 & & 0.788 & 1\\
& $d=100$ & 0.855 & 1.000 & & 0.904 & 1.000 & & 0.653 & 1\\\hline
EE$_{\Omega}$ & $d=50$ & 0.510 & 0.436 & & 0.503 & 0.435 & & 0.451 & 0.407\\
& $d=100$ & 0.481 & 0.421 & & 0.473 & 0.419 & & 0.433 & 0.397\\\hline
\end{tabular}
\end{table}
\medskip
\noindent\textbf{Example 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{tab:3} 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{tab:4} 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{tab:3} 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$.
\begin{table}
\caption{Transition matrix and Granger network estimation in Example 2.}
\label{tab:3}
\centering
\begin{tabular}
[c]{llllllllll}\hline\hline
& & \multicolumn{2}{l}{tv-wgLASSO} & & \multicolumn{2}{l}{tv-Oracle} & &
\multicolumn{2}{l}{tv-Full}\\\cline{3-4}\cline{6-7}\cline{9-10}
measure & dimension & $n=200$ & $n=400$ & & $n=200$ & $n=400$ & & $n=200$ &
$n=400$\\\hline
FP & $d=50$ & 13.53 & 12.75 & & 0 & 0 & & 2401 & 2401\\
& $d=100$ & 20.73 & 7.73 & & 0 & 0 & & - & 9801\\
FN & $d=50$ & 18.56 & 11.11 & & 0 & 0 & & 0 & 0\\
& $d=100$ & 37.55 & 13.90 & & 0 & 0 & & - & 0\\
TPR & $d=50$ & 0.813 & 0.888 & & 1 & 1 & & 1 & 1\\
& $d=100$ & 0.811 & 0.930 & & 1 & 1 & & - & 1\\
TNR & $d=50$ & 0.994 & 0.995 & & 1 & 1 & & 0 & 0\\
& $d=100$ & 0.998 & 0.999 & & 1 & 1 & & - & 0\\
PPV & $d=50$ & 0.859 & 0.875 & & 1 & 1 & & 0.040 & 0.040\\
& $d=100$ & 0.888 & 0.960 & & 1 & 1 & & - & 0.020\\
NPV & $d=50$ & 0.992 & 0.995 & & 1 & 1 & & 0 & 0\\
& $d=100$ & 0.996 & 0.999 & & 1 & 1 & & - & 0\\
F1 & $d=50$ & 0.834 & 0.881 & & 1 & 1 & & 0.076 & 0.076\\
& $d=100$ & 0.847 & 0.945 & & 1 & 1 & & - & 0.039\\
MCC & $d=50$ & 0.828 & 0.876 & & 1 & 1 & & 0 & 0\\
& $d=100$ & 0.846 & 0.943 & & 1 & 1 & & - & 0\\\hline
average $R^{2}$ & $d=50$ & 0.465 & 0.448 & & 0.477 & 0.462 & & 0.963 &
0.829\\
& $d=100$ & 0.473 & 0.467 & & 0.483 & 0.471 & & - & 0.978\\\hline
EE$_{A}$ & $d=50$ & 0.328 & 0.250 & & 0.171 & 0.122 & & 58.44 & 1.510\\
& $d=100$ & 0.323 & 0.204 & & 0.168 & 0.122 & & - & 82.60\\
RMSE$_{e}$ & $d=50$ & 0.631 & 0.476 & & 0.417 & 0.305 & & 1.673 & 1.414\\
& $d=100$ & 0.613 & 0.390 & & 0.414 & 0.309 & & - & 1.720\\\hline
\end{tabular}
\end{table}
\begin{table}
\caption{Precision matrix and partial correlation network estimation in
Example 2.}
\label{tab:4}
\centering
\begin{tabular}
[c]{llllllllll}\hline\hline
& & \multicolumn{2}{l}{tv-CLIME} & & \multicolumn{2}{l}{infeasible tv-CLIME}
& & \multicolumn{2}{l}{tv-GLASSO}\\\cline{3-4}\cline{6-7}\cline{9-10}
measure & dimension & $n=200$ & $n=400$ & & $n=200$ & $n=400$ & & $n=200$ &
$n=400$\\\hline
FP & $d=50$ & 0.03 & 0.04 & & 0.02 & 0.03 & & 0 & 0.01\\
& $d=100$ & 0.01 & 0 & & 0 & 0.01 & & 0 & 0.01\\
FN & $d=50$ & 12.62 & 0.82 & & 2.34 & 0 & & 20.84 & 0.06\\
& $d=100$ & 24.71 & 0.23 & & 6.21 & 0.01 & & 49.73 & 0.43\\
TPR & $d=50$ & 0.742 & 0.983 & & 0.952 & 1 & & 0.575 & 0.997\\
& $d=100$ & 0.750 & 0.998 & & 0.937 & 1.000 & & 0.498 & 0.996\\
TNR & $d=50$ & 1.000 & 1.000 & & 1.000 & 1.000 & & 1 & 1.000\\
& $d=100$ & 1.000 & 1 & & 1 & 1.000 & & 1 & 1.000\\
PPV & $d=50$ & 0.999 & 0.999 & & 1.000 & 0.999 & & 1 & 1.000\\
& $d=100$ & 1.000 & 1 & & 1 & 1.000 & & 1 & 1.000\\
NPV & $d=50$ & 0.989 & 0.999 & & 0.998 & 1 & & 0.983 & 1.000\\
& $d=100$ & 0.995 & 1.000 & & 0.999 & 1.000 & & 0.990 & 1.000\\
F1 & $d=50$ & 0.850 & 0.991 & & 0.975 & 1.000 & & 0.725 & 0.998\\
& $d=100$ & 0.857 & 0.999 & & 0.967 & 1.000 & & 0.662 & 0.998\\
MCC & $d=50$ & 0.856 & 0.991 & & 0.975 & 1.000 & & 0.749 & 0.998\\
& $d=100$ & 0.864 & 0.999 & & 0.967 & 1.000 & & 0.701 & 0.998\\\hline
EE$_{\Omega}$ & $d=50$ & 0.598 & 0.533 & & 0.526 & 0.485 & & 0.560 & 0.514\\
& $d=100$ & 0.560 & 0.489 & & 0.486 & 0.458 & & 0.536 & 0.496\\\hline
\end{tabular}
\end{table}
\medskip
\noindent\textbf{Example 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{tab:5} 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.
\begin{table}
\caption{Estimation accuracy of dual networks in Example 3.}
\label{tab:5}
\centering
\begin{tabular}
[c]{llllllllll}\hline\hline
& & \multicolumn{2}{c}{tv-wgLASSO} & & \multicolumn{2}{c}{tv-Oracle} & &
\multicolumn{2}{c}{tv-Full}\\\cline{3-4}\cline{6-7}\cline{9-10}
measure & dimension & $n=200$ & $n=400$ & & $n=200$ & $n=400$ & & $n=200$ &
$n=400$\\\hline
average $R^{2}$ & $d=50$ & 0.009 & 0.029 & & 0.891 & 0.588 & & 0.891 &
0.588\\
& $d=100$ & 0.005 & 0.020 & & - & 0.930 & & - & 0.930\\
EE$_{A}$ & $d=50$ & 0.383 & 0.348 & & 56.66 & 1.927 & & 56.66 & 1.927\\
& $d=100$ & 0.388 & 0.364 & & - & 97.60 & & - & 97.60\\
RMSE$_{e}$ & $d=50$ & 0.515 & 0.463 & & 1.716 & 1.300 & & 1.716 & 1.300\\
& $d=100$ & 0.523 & 0.486 & & - & 1.776 & & - & 1.776\\\hline
& & \multicolumn{2}{c}{tv-CLIME} & & \multicolumn{2}{c}{infeasible tv-CLIME}
& & \multicolumn{2}{c}{tv-GLASSO}\\\cline{3-4}\cline{6-7}\cline{9-10}
& & $n=200$ & $n=400$ & & $n=200$ & $n=400$ & & $n=200$ & $n=400$\\\hline
EE$_{\Omega}$ & $d=50$ & 1.669 & 1.601 & & 1.613 & 1.572 & & 1.584 & 1.570\\
& $d=100$ & 1.674 & 1.615 & & 1.616 & 1.580 & & 1.587 & 1.588\\\hline
\end{tabular}
\end{table}
\medskip
\noindent\textbf{Example 4}.\ \ The data is generated from a factor-adjusted
time-varying VAR model in the form of (\ref{eq5.2}). 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{tab:6} reports the estimation results of the time-varying
transition matrices and Granger networks for the idiosyncratic errors, and
Table \ref{tab:7} reports the estimation results of the time-varying precision
matrices and partial correlation networks. Comparing with the results in
Tables \ref{tab:3} and \ref{tab:4}, 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$).
\begin{table}
\begin{minipage}{0.45\linewidth}
\caption{\label{tab:6}Factor-adjusted transition matrix and Granger network estimation in Example 4.}\centering
\begin{tabular}{llll}\hline\hline
&& \multicolumn{2}{l}{tv-wgLASSO} \\ \cline{3-4}
measure&dimension&$n=200$ &$n=400$ \\ \hline
FP &$d=50$&11.35&10.60 \\
& $d=100$&20.40&10.41 \\
FN & $d=50$&35.97&14.77 \\
& $d=100$& 65.45&20.68\\
TPR& $d=50$&0.637&0.851\\
& $d=100$&0.671 &0.896 \\
TNR& $d=50$&0.995&0.996 \\
& $d=100$&0.998&0.999 \\
PPV& $d=50$&0.852&0.890 \\
& $d=100$& 0.869 &0.945 \\
NPV& $d=50$&0.985&0.994 \\
& $d=100$&0.993&0.998 \\
F1& $d=50$&0.725&0.869 \\
& $d=100$&0.756&0.920 \\
MCC& $d=50$&0.725&0.865 \\
& $d=100$&0.759&0.919 \\\hline
average $R^2$&$d=50$&0.298&0.350 \\
& $d=100$& 0.339&0.389 \\\hline
EE$_A$&$d=50$& 0.413&0.283 \\
& $d=100$&0.396&0.241 \\
RMSE$_e$ &$d=50$&1.319&1.025 \\
& $d=100$&1.230&0.856 \\\hline
\end{tabular}
\end{minipage}
\hfill\begin{minipage}{0.45\linewidth}
\caption{\label{tab:7}Factor-adjusted precision matrix and partial correlation network estimation in Example 4.}\centering
\begin{tabular}{llll}\hline\hline
&& \multicolumn{2}{l}{tv-CLIME} \\ \cline{3-4}
measure&dimension&$n=200$ &$n=400$ \\ \hline
FP &$d=50$&0.01&0.01 \\
& $d=100$&0&0.02\\
FN & $d=50$&38.22&5.36 \\
& $d=100$&65.99&2.21 \\
TPR& $d=50$&0.220&0.891\\
& $d=100$&0.333&0.978 \\
TNR& $d=50$&1.000&1.000 \\
& $d=100$&1&1.000 \\
PPV& $d=50$&0.999&1.000 \\
& $d=100$&1&1.000 \\
NPV& $d=50$&0.969&0.995 \\
& $d=100$&0.987&1.000 \\
F1& $d=50$&0.349&0.941 \\
& $d=100$&0.496& 0.989 \\
MCC& $d=50$&0.448&0.941 \\
& $d=100$&0.570&0.988 \\\hline
EE$_\Omega$&$d=50$&0.670&0.585 \\
& $d=100$&0.628&0.534 \\
\hline
\end{tabular}
\end{minipage}
\end{table}
\section{An empirical application}
\label{sec7} \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 \cite{MN16}.
We follow \cite{MN16} and \cite{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{sec5} to accommodate
strong cross-sectional dependence: the approximate factor model (\ref{eq5.1})
with constant factor loadings, and the time-varying factor model (\ref{eq5.2})
with dynamic factor loadings. The information criteria proposed by \cite{BN02}
and \cite{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.
\cite{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}}$.
\begin{figure}
\begin{center}
\includegraphics[scale=0.3]{pic/VAR1mgmCnetDynfac.png}
\includegraphics[scale=0.3]{pic/VAR1netCnetDynfac.png}
\includegraphics[scale=0.3]{pic/VAR1mgmCmatDynfac.png}
\includegraphics[scale=0.3]{pic/VAR1netCmatDynfac.png}
\end{center}
\caption{{\protect\small The estimated Granger causality networks using the
factor-adjusted static VAR(1) model (left) and time-varying VAR(1) model
(right). }}
\label{fig2}
\end{figure}
Figure \ref{fig2} 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{fig3} 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 \cite{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.
\begin{figure}
\begin{center}
\includegraphics[scale=0.3]{pic/VAR1mgmCnetNofac.png}
\includegraphics[scale=0.3]{pic/VAR1netCnetNofac.png}
\includegraphics[scale=0.3]{pic/VAR1mgmCmatNofac.png}
\includegraphics[scale=0.3]{pic/VAR1netCmatNofac.png}
\end{center}
\caption{{\protect\small The estimated Granger causality networks using the
static VAR(1) model (left) and time-varying VAR(1) model (right) without
factor-adjustment. }}
\label{fig3}
\end{figure}
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 ``\textsf{tvReg}". Figure \ref{fig4} 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.
\begin{figure}[tbh]
\begin{center}
\includegraphics[scale=0.35]{pic/CoeffBOGMBASE.png}
\includegraphics[scale=0.35]{pic/CoeffSP500.png}
\includegraphics[scale=0.35]{pic/CoeffSP500ind.png}
\includegraphics[scale=0.35]{pic/CoeffSPPE.png}
\end{center}
\caption{{\protect\small The estimated time-varying coefficients linked to
S\&P PE ratio with 90\% confident bands.}}
\label{fig4}
\end{figure}
We plot the estimated partial correlation networks in Figure \ref{fig5}, 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 \citep{JV15}, and draw the 90\% confidence bands using
R package ``\textsf{SILGGM}" in Figure \ref{fig6}. 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{fig6} suggest that
time-invariant partial correlation linkages are inappropriate to describe the
network structure of the FRED-MD data.
\begin{figure}[tbh]
\begin{center}
\includegraphics[scale=0.3]{pic/VAR1PnetDynfac.png}
\includegraphics[scale=0.3]{pic/VAR1PnetNofac.png}
\includegraphics[scale=0.3]{pic/VAR1PmatDynfac.png}
\includegraphics[scale=0.3]{pic/VAR1PmatNofac.png}
\end{center}
\caption{{\protect\small The estimated partial correlation networks with
(left) and without (right) factor adjustment. }}
\label{fig5}
\end{figure}
\smallskip
\begin{figure}[tbh]
\begin{center}
\includegraphics[scale=0.35]{pic/PnetCoeffSP500.png}
\includegraphics[scale=0.35]{pic/PnetCoeffSP500ind.png}
\includegraphics[scale=0.35]{pic/PnetCoeffSPdiv.png}
\includegraphics[scale=0.35]{pic/PnetCoeffBAAFFM.png}
\end{center}
\caption{{\protect\small The estimated time-varying elements in the precision
matrix linked to S\&P PE ratio with 90\% confident bands.}}
\label{fig6}
\end{figure}
\section{Conclusion}
\label{sec8}
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.
\section*{\Large Supplementary materials}
{\small The supplement contains proofs of the main asymptotic theorems, some
technical lemmas with proofs, verification of Assumption \ref{ass:3}(ii) and
discussions on tuning parameter selection.}
\begin{thebibliography}{}
{\footnotesize \harvarditem{Bai \harvardand\ Ng}{2002}{BN02} \textsc{Bai, J.
\harvardand\ Ng, S.} (2002). Determining the number of factors in approximate
factor models. \emph{Econometrica} 90, 191--221. }
{\footnotesize \harvarditem{Barigozzi \harvardand\ Brownlees}{2019}{BB19}
\textsc{Barigozzi, M. \harvardand\ Brownlees, C.} (2019). NETS: Network
estimation for time series. \emph{Journal of Applied Econometrics} 34,
347--364. }
{\footnotesize \harvarditem{Barigozzi, Cho \harvardand\ Owens}{2022}{BCO22}
\textsc{Barigozzi, M., Cho, H. and Owens, D.} (2022). FNETS: Factor-adjusted
network estimation and forecasting for high-dimensional time series. Working
paper available at \url{https://arxiv.org/pdf/2201.06110.pdf}. }
{\footnotesize \harvarditem{Basu \harvardand\ Michailidis}{2015}{BM15}
\textsc{Basu, S. \harvardand\ Michailidis, G.} (2015). Regularized estimation
in sparse high-dimensional time series models. \emph{The Annals of Statistics}
43, 1535--1567. }
{\footnotesize \harvarditem{Basu, Shojaie \harvardand\ Michailidis}{2015}{BSM15}
\textsc{Basu, S. Shojaie, A. \harvardand\ Michailidis, G.} (2015). Network
Granger causality with inherent grouping structure. \emph{Journal of Machine
Learning Research} 16, 417--453. }
{\footnotesize \harvarditem{Bickel, Ritov \harvardand\ Tsybakov}{2009}{BRT09}
\textsc{Bickel, P., Ritov, Y. and Tsybakov, A.} (2009). Simultaneous analysis
of lasso and dantzig selector. \emph{The Annals of Statistics}, 37,
1705--1732. }
{\footnotesize \harvarditem{Burt, Kilduff \harvardand\ Tasselli}{2013}{BKT13}
\textsc{Burt, R. S., Kilduff, M. and Tasselli, S.} (2013). Social network
analysis: foundations and frontiers on advantage. \emph{Annual Review of
Psychology} 64, 527--547. }
{\footnotesize
}
{\footnotesize \harvarditem{Cai, Liu \harvardand\ Luo}{2011}{CLL11}
\textsc{Cai, T. T., Liu, W. \harvardand\ Luo, X.} (2011). A constrained
$\ell_{1}$ minimization approach to sparse precision matrix estimation.
\emph{Journal of the American Statistical Association} 106, 594--607. }
{\footnotesize \harvarditem{Cai}{2007}{C07} \textsc{Cai, Z.} (2007). Trending
time-varying coefficient time series models with serially correlated errors.
\emph{Journal of Econometrics} 136, 163--188. }
{\footnotesize \harvarditem{Chamberlain \harvardand\ Rothschild}{1983}{CR83}
\textsc{Chamberlain, G. \harvardand\ Rothschild, M.} (1983). Arbitrage, factor
structure and mean-variance analysis in large asset markets.
\emph{Econometrica} 51, 1305--1324. }
{\footnotesize \harvarditem{Chen, Fan \harvardand\ Zhu}{2020}{CFZ20}
\textsc{Chen, E., Fan, J. \harvardand\ Zhu, X.} (2020). Community network
autoregression for high-dimensional time series. Working paper available at
\url{https://arxiv.org/abs/2007.05521}. }
{\footnotesize \harvarditem{Chen {\em et al}.}{2021}{CLWZ21} \textsc{Chen, J.,
Li, D., Wei, L. and Zhang, W.} (2021). Nonparametric homogeneity pursuit in
functional-coefficient models. \emph{Journal of Nonparametric Statistics} 33,
387--416. }
{\footnotesize \harvarditem{Cheng {\em et al}.}{2014}{CHLP14} \textsc{Cheng,
M., Honda, T., Li, J. \harvardand\ Peng, H.} (2014). Nonparametric
independence screening and structure identification for ultra-high dimensional
longitudinal data. \emph{The Annals of Statistics}, 42, 1819--1849. }
{\footnotesize \harvarditem{Dahlhaus}{1997}{D97} \textsc{Dahlhaus, R.} (1997).
Fitting time series models to nonstationary processes. \emph{The Annals of
Statistics} 25, 1--37. }
{\footnotesize \harvarditem{Dahlhaus \harvardand\ Subba Rao}{2006}{DS06}
\textsc{Dahlhaus, R. \harvardand\ Subba Rao, S.} (2006). Statistical inference
for time-varying ARCH processes. \emph{The Annals of Statistics} 34,
1075--1114. }
{\footnotesize \harvarditem{Davis, Zang \harvardand\ Zheng}{2016}{DZZ16}
\textsc{Davis, R., Zang, P. \harvardand\ Zheng, T.} (2016). Sparse vector
autoregressive modeling. \emph{Journal of Computational and Graphical
Statistics} 25, 1077--1096. }
{\footnotesize \harvarditem{Dempster}{1972}{D72} \textsc{Dempster, A.P.}
(1972). Covariance selection. \emph{Biometrics} 28, 157--175. }
{\footnotesize \harvarditem{Diebold \harvardand\ Ylmaz}{2014}{DY14}
\textsc{Diebold, F. \harvardand\ Yilmaz, K.} (2014). On the network topology
of variance decompositions: Measuring the connectedness of financial firms.
\emph{Journal of Econometrics} 182, 119--134. }
{\footnotesize \harvarditem{Diebold \harvardand\ Ylmaz}{2015}{DY15}
\textsc{Diebold, F. \harvardand\ Yilmaz, K.} (2015). \emph{Financial and
Macroeconomic Connectedness: A Network Approach to Measurement and
Monitoring}. Oxford University Press. }
{\footnotesize \harvarditem{Ding, Qiu \harvardand\ Chen}{2017}{DQC17}
\textsc{Ding, X., Qiu, Z. \harvardand\ Chen, X.} (2017). Sparse transition
matrix estimation for high-dimensional and locally stationary vector
autoregressive models. \emph{Electronic Journal of Statistics} 11, 3871--3902.
}
{\footnotesize \harvarditem{Fan, Feng \harvardand\ Wu}{2009}{FFW09}
\textsc{Fan, J., Feng, Y. \harvardand\ Wu, Y.} (2009). Network exploration via
the adaptive lasso and SCAD penalties. \emph{The Annals of Applied Statistics}
3, 521--541. }
{\footnotesize \harvarditem{Fan \harvardand\ Gijbels}{1996}{FG96} \textsc{Fan,
J. and Gijbels, I.} (1996). \emph{Local Polynomial Modelling and Its
Applications}. Chapman \& Hall. }
{\footnotesize \harvarditem{Fan, Masini \harvardand\ Medeiros}{2021}{FMM21}
\textsc{Fan, J., Masini, R. and Medeiros, M.} (2021). Bridging factor and
sparse models. Working paper available at
\url{https://arxiv.org/abs/2102.11341}. }
{\footnotesize \harvarditem{Fan \harvardand\ Li}{2001}{FL01} \textsc{Fan, J.
\harvardand\ Li, R.} (2001). Variable selection via nonconcave penalized
likelihood and its oracle properties. \emph{Journal of the American
Statistical Association} 96, 1348--1360. }
{\footnotesize \harvarditem{Fan, Liao \harvardand\ Mincheva}{2013}{FLM13}
\textsc{Fan, J., Liao, Y. and Mincheva, M.} (2013). Large covariance
estimation by thresholding principal orthogonal complements (with discussion).
\emph{Journal of the Royal Statistical Society, Series B} 75, 603--680. }
{\footnotesize
}
{\footnotesize \harvarditem{Fan, Ma \harvardand\ Dai}{2014}{FMD14}
\textsc{Fan, J., Ma, Y. \harvardand\ Dai, W.} (2014). Nonparametric
independence screening in sparse ultra-high dimensional varying coefficient
models. \emph{Journal of the American Statistical Association} 109,
1270--1284. }
{\footnotesize \harvarditem{Fan, Xue \harvardand\ Zou}{2014}{FXZ14}
\textsc{Fan, J., Xue, L. \harvardand\ Zou, H.} (2014). Strong oracle
optimality of folded concave penalized estimation. \emph{The Annals of
Statistics} 42, 819--849. }
{\footnotesize \harvarditem{Granger}{1969}{G69} \textsc{Granger, C. W.}
(1969). Investigating causal relations by econometric models and
cross-spectral methods. \emph{Econometrica} 37, 424--438. }
{\footnotesize \harvarditem{Hafner \harvardand\ Linton}{2010}{HL10}
\textsc{Hafner, C. \harvardand\ Linton, O.} (2010). Efficient estimation of a
multivariate multiplicative volatility model. \emph{Journal of Econometrics}
159, 55--73. }
{\footnotesize \harvarditem{Han, Lu \harvardand\ Liu}{2015}{HLL15}
\textsc{Han, F., Lu, H. \harvardand\ Liu H.} (2015). A direct estimation of
high dimensional stationary vector autoregressions. \emph{Journal of Machine
Learning Research} 16, 3115--3150. }
{\footnotesize \harvarditem{Hautsch, Schaumburg \harvardand\ Schienle}{2014}{HSS14}
\textsc{Hautsch, N., Schaumburg, J. \harvardand\ Schienle, M.} (2014).
Forecasting systemic impact in financial networks. \emph{International Journal
of Forecasting} 30, 781--794. }
{\footnotesize \harvarditem{Jankova \harvardand\ van de Geer}{2015}{JV15}
\textsc{Jankova, J. \harvardand\ van de Geer S.} (2015). Confidence intervals
for high-dimensional inverse covariance estimation. \emph{Electronic Journal
of Statistics} 9, 1205--1229. }
{\footnotesize
}
{\footnotesize \harvarditem{Kock \harvardand\ Callot}{2015}{KC15}
\textsc{Kock, A.B. and Callot, L.} (2015). Oracle inequalities for high
dimensional vector autoregressions. \emph{Journal of Econometrics} 186,
325--344. }
{\footnotesize \harvarditem{Kolar {\em et al}.}{2010}{KSAX10} \textsc{Kolar,
M., Song, L. Ahmed, A. and Xing, E.} (2010). Estimating time-varying networks.
\emph{The Annals of Applied Statistics} 4, 94--123. }
{\footnotesize \harvarditem{Koo \harvardand\ Linton}{2012}{KL12} \textsc{Koo,
B. and Linton O.} (2012). Estimation of semiparametric locally stationary
diffusion models. \emph{Journal of Econometrics} 170, 210--233. }
{\footnotesize \harvarditem{Krampe \harvardand\ Margaritella}{2022}{KM22}
Krampe, J. and Margaritella, L. (2022). Factor models with sparse VAR
idiosyncratic components. Working paper available at
\url{https://arxiv.org/pdf/2112.07149.pdf}. }
{\footnotesize \harvarditem{Lam \harvardand\ Fan}{2009}{LF09} \textsc{Lam, C.
\harvardand\ Fan, J.} (2009). Sparsity and rates of convergence in large
covariance matrix estimation. \emph{The Annals of Statistics} 37, 4254--4278.
}
{\footnotesize \harvarditem{Li, Chen \harvardand\ Gao}{2011}{LCG11}
\textsc{Li, D., Chen, J. \harvardand\ Gao, J.} (2011). Nonparametric
time-varying coefficient panel data models with fixed effects. \emph{The
Econometrics Journal} 14, 387--408. }
{\footnotesize \harvarditem{Li, Ke and Zhang}{2015}{LKZ15} \textsc{Li, D., Ke,
Y. \harvardand\ Zhang, W.} (2015). Model selection and structure specification
in ultra-high dimensional generalised semi-varying coefficient models.
\emph{The Annals of Statistics} 43, 2676--2705. }
{\footnotesize
}
{\footnotesize \harvarditem{Lian}{2012}{L12} \textsc{Lian, H.} (2012).
Variable selection for high-dimensional generalized varying-coefficient
models. \emph{Statistica Sinica}, 22, 1563--1588. }
{\footnotesize \harvarditem{Liu, Li and Wu}{2014}{LLW14} \textsc{Liu, J., Li,
R. and Wu, R.} (2014). Feature selection for varying coefficient models with
ultrahigh dimensional covariates. \emph{Journal of the American Statistical
Association}, 109, 266--274. }
{\footnotesize \harvarditem{Loh \harvardand\ Wainwright}{2013}{LW13}
\textsc{Loh, P. \harvardand\ Wainwright, M.} (2013). Structural estimation for
discrete graphical models: Generalized covariance matrices and their inverse.
\emph{The Annals of Statistics} 41, 3022--3049. }
{\footnotesize \harvarditem{Liu \harvardand\ Zhang}{2021}{LZ21}\textsc{Liu, L.
\harvardand\ Zhang, D.} (2021). Robust estimation of high-dimensional vector
autoregressive models. Working paper available at
\url{https://arxiv.org/abs/2109.10354}. }
{\footnotesize \harvarditem{L\"utkepohl}{2006}{Lu06} \textsc{L\"utkepohl, H.}
(2006). \emph{New Introduction to Multiple Time Series Analysis}. Springer. }
{\footnotesize \harvarditem{McCracken \harvardand\ Ng}{2016}{MN16}
\textsc{McCracken, M.W. \harvardand\ Ng, S.} (2016). FRED-MD: A monthly
database for macroeconomic research. \emph{Journal of Business \& Economic
Statistics} 34, 574--589. }
{\footnotesize \harvarditem{McCracken \harvardand\ Ng}{2020}{MN20}
\textsc{McCracken, M.W. \harvardand\ Ng, S.} (2020). FRED-QD: A quarterly
database for macroeconomic research. Working paper available at
\url{https://www.nber.org/papers/w26872}. }
{\footnotesize \harvarditem{Miao, Phillips \harvardand\ Su}{2022}{MPS22}
\textsc{Miao, K., Phillips, P.C.B. \harvardand\ Su, L.} (2022).
High-dimensional VARs with common factors. Forthcoming in \emph{Journal of
Econometrics}. }
{\footnotesize \harvarditem{Motta, Hafner \harvardand\ von Sachs}{2011}{MHvS11}
Motta, G., Hafner, C. and von Sachs, R. (2011). Locally stationary factor
models: identification and nonparametric estimation. \emph{Econometric Theory}
27, 1279--1319. }
{\footnotesize \harvarditem{Newman}{2002}{N02} \textsc{Newman, M. E. J.}
(2002). Spread of epidemic disease on networks. \emph{Physics Review, Series
E} 66, 016128. }
{\footnotesize \harvarditem{Safikhani \harvardand\ Shojaie}{2022}{SS22}
\textsc{Safikhani, A. \harvardand\ Shojaie, A.} (2022). Joint structural break
detection and parameter estimation in high-dimensional non-stationary VAR
models. \emph{Journal of the American Statistical Association} 117, 251--264.
}
{\footnotesize \harvarditem{Scott}{2017}{S17} \textsc{Scott, J.} (2017).
\emph{Social Network Analysis} (4th Edition). Sage, London. }
{\footnotesize \harvarditem{Stock \harvardand\ Watson}{2002}{SW02} Stock, J.
H. and Watson, M. W. (2002). Forecasting using principal components from a
large number of predictors. \emph{Journal of the American Statistical
Association} 97, 1167--1179. }
{\footnotesize \harvarditem{Su \harvardand\ Wang}{2017}{SW17} \textsc{Su, L.
and Wang, X.} (2017). On time-varying factor models: estimation and testing.
\emph{Journal of Econometrics} 198, 84--101. }
{\footnotesize
}
{\footnotesize
}
{\footnotesize \harvarditem{Tibshirani}{1996}{T96} \textsc{Tibshirani, R. J.}
(1996). Regression shrinkage and selection via the LASSO. \emph{Journal of the
Royal Statistical Society Series B} 58, 267--288. }
{\footnotesize \harvarditem{Vogt}{2012}{Vo12} \textsc{Vogt, M.} (2012).
Nonparametric regression for locally stationary time series. \emph{The Annals
of Statistics} 40, 2601--2633. }
{\footnotesize \harvarditem{Wainwright}{2019}{W19} \textsc{Wainwright, M. J.}
(2019). \emph{High-Dimensional Statistics: A Non-Asymptotic Viewpoint}.
Cambridge Series in Statistical and Probabilistic Mathematics. }
{\footnotesize \harvarditem{Wang, Li and Huang}{2008}{WLH08} \textsc{Wang, L.,
Li, H. and Huang, J.} (2008). Variable selection in nonparametric
varying-coefficient models for analysis of repeated measurements.
\emph{Journal of the American Statistical Association} 103, 1556--1569. }
{\footnotesize \harvarditem{Wang and Xia}{2009}{WX09} \textsc{Wang, H. and
Xia, Y.} (2009). Shrinkage estimation of the varying-coefficient model.
\emph{Journal of the American Statistical Association} 104, 747--757. }
{\footnotesize \harvarditem{Wang, Yu and Rinaldo}{2021}{WYR21} \textsc{Wang,
D., Yu, Y. \harvardand\ Rinaldo, A.} (2021). Optimal change point detection
and localization in sparse dynamic networks. \emph{The Annals of Statistics}
49, 203--232. }
{\footnotesize \harvarditem{Xu, Chen \harvardand\ Wu}{2020}{XCW20} \textsc{Xu,
M., Chen, X. and Wu, W.} (2020). Estimation of dynamic networks for
high-dimensional nonstationary time series. \emph{Entropy} 22, 55. }
{\footnotesize \harvarditem{Yan, Gao \harvardand\ Peng}{2020}{YGP20}
\textsc{Yan, Y., Gao, J. \harvardand\ Peng, B.} (2020). A class of
time-varying vector moving average $(\infty)$ models. Working paper available
at \url{https://arxiv.org/abs/2010.01492}. }
{\footnotesize \harvarditem{Yuan}{2010}{Y10} \textsc{Yuan, M.} (2010). High
dimensional inverse covariance matrix estimation via linear programming.
\emph{Journal of Machine Learning Research} 11, 2261--2286. }
{\footnotesize \harvarditem{Yuan \harvardand\ Lin}{2007}{YL07} \textsc{Yuan,
M. \harvardand\ Lin, Y.} (2007). Model selection and estimation in the
Gaussian graphical model. \emph{Biometrika} 94, 19--35. }
{\footnotesize \harvarditem {Zhang \harvardand\ Wu}{2012}{ZW12} \textsc{Zhang,
T. and Wu, W. B.} (2012). Inference of time varying regression models.
\emph{The Annals of Statistics} 40, 1376--1402. }
{\footnotesize \harvarditem {Zhang \harvardand\ Wu}{2021}{ZW21} \textsc{Zhang,
D. and Wu, W.} (2021). Convergence of covariance and spectral density
estimators for high-dimensional locally stationary processes. \emph{The Annals
of Statistics} 49, 233--254. }
{\footnotesize \harvarditem{Zhao {\em et al}}{2022}{ZLWL22} \textsc{Zhao, J.,
Liu, X., Wang, H. and Leng, C.} (2022). Dimension reduction for covariates in
network data. \emph{Biometrika} 109, 85--102. }
{\footnotesize
}
{\footnotesize \harvarditem{Zhou, Lafferty \harvardand\ Wasserman}{2010}{ZLW10}
\textsc{Zhou, S., Lafferty, J. and Wasserman, L.} (2010). Time varying
undirected graphs. \emph{Machine Learning} 80, 295--319. }
{\footnotesize \harvarditem{Zhu {\em et al}.}{2019}{ZCLW19} \textsc{Zhu, X.,
Chang, X., Li, R. \harvardand\ Wang, H.} (2019). Portal nodes screening for
large scale social networks. \emph{Journal of Econometrics} 209, 145--157. }
{\footnotesize \harvarditem{Zhu {\em et al}.}{2017}{ZPLLW17} \textsc{Zhu, X.,
Pan, R., Li, G., Liu, Y. \harvardand\ Wang, H.} (2017). Network vector
autoregression. \emph{The Annals of Statistics} 45, 1096--1123. }
{\footnotesize \harvarditem{Zou \harvardand\ Li}{2008}{ZL08} \textsc{Zou, H.
and Li, R.} (2008). One-step sparse estimates in nonconcave penalized
likelihood models (with discussion). \emph{The Annals of Statistics}, 36,
1509--1566. }
\end{thebibliography}
\newpage
\begin{center}
{\LARGE\bf Supplement to ``Estimating Time-Varying Networks for High-Dimensional Time Series"}
\end{center}
\maketitle