Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
56,560 characters · 20 sections · 50 citation commands
Learning Time-Varying Correlation Networks with FDR Control via Time-Varying P-values
Keywords: time-varying correlation network, non-stationary time series, Gaussian approximation, multiplier bootstrap, false discovery rate
Recent years have seen rapidly increasing interest in time-varying brain functional connectivity lurie2020questions, jiang2022time. Typically, brain signals are represented as time series, and the most widely used measure of functional connectivity is the correlation between these time series lurie2020questions. Similarly, in fields such as risk management DIEBOLD2014119, Marti2021, correlation relationships between variables often evolve due to changing underlying processes, external influences, or latent factors. In these contexts, time-varying correlation networks provide a powerful framework to analyze how interactions between variables change over time.
Figure (ref) illustrates a simple example of constructing time-varying correlation networks from an electroencephalogram (EEG) data set in zhang1997electrophysiological used in our real data analysis. In these networks, two vertices are connected whenever the null hypothesis of no pairwise correlation is rejected at a given time point. Existing approaches for constructing (time-varying) correlation networks frequently neglect proper uncertainty control; see, for example, kazemilari2015correlation,bailey2019multiple,masuda2025introduction. Effective uncertainty control is especially important in exploratory research, where the goal is to maximize the number of genuine discoveries while remaining robust to the false positive paradox, in which a considerable proportion of reported findings may in fact be spurious.
A widely adopted criterion for controlling uncertainty is the False Discovery Rate (FDR), which aims to limit the expected proportion of false positives (incorrect rejections of the null hypothesis) among all discoveries (rejected null hypotheses), as discussed in benjamini2001control. An alternative criterion is Family-Wise Error Rate (FWER), which is the probability of making any false discoveries. In contrast to the FWER control which can lead to overly conservative inferences in large scale data, FDR control offers a more balanced approach and is thus extensively utilized across various scientific disciplines.
Given that the time-varying correlation network is defined at infinitely many time points, our goal is to control the FDR uniformly in the sense that \[ \lim_{n\rightarrow\infty}\sup_{t\in[b_n,1-b_n]} \left\{ \mathrm{FDR}\left( t \right) - \alpha \frac{2|{\cal H}_0(t)|}{p(p-1)} \right\} \le 0, \] where $\alpha$ is the prescribed FDR control level, $\mathrm{FDR}\left( t \right)$ denotes the FDR of the correlation network at time $t$, $2|{\cal H}_0(t)|/p(p-1)$ is the proportion of true nulls at time $t$, and $b_n$ is a diminishing sequence. The supremum over $t$ represents uniform FDR control across uncountably many time points over the whole interval, which adds a significant challenge to our framework.
Moreover, the increasing complexity of data poses significant challenges to existing analytical frameworks. As shown in Figure (ref), the data exhibit potential mean discontinuities, non-stationarity, and high dimensionality, which are also traits summarized in lurie2020questions and jiang2022time. Additionally, non-Gaussianity is especially relevant in financial data analysis, see de2023modeling and jondeau2007financial. These challenges highlight the need for a new robust method and a theoretical framework in time-varying correlation network construction.
Outline and contributions. In this paper, we propose a new framework to learn time-varying correlation networks with uniformly controlled FDR from locally stationary non-Gaussian high-dimensional time series with complex trends. Our major contributions are:
Related work. The problem of multiple testing of correlations in large-scale data with uncertainty control has been investiated by cai2016large and bailey2019multiple. Recently, learning time-varying network structures from high-dimensional data has gained increasing attention. basu2023graphical and CHEN2025network developed advanced methods for high-dimensional data width dependence and changing temporal dynamics. However, the former can only construct a static network from the time series data, while the latter focuses on inferring time-varying networks from partial correlations and without providing guaranties for uncertainty control. These works also ignore potential jumps in the mean of observations. Meanwhile, bai2025time proposed a framework for learning time-varying correlation network structures, but via FWER control instead of FDR control.
The foundational work on FDR control was established by benjamini1995controlling and benjamini2001control, who introduced the Benjamini--Hochberg (B-H) procedure with independent P-values and the Benjamini--Yekutieli (B-Y) procedure with arbitrarily dependent P-values for FDR control. There are other studies addressing dependent P-values, such as barber2015controlling and fithian2022conditional. However, the complex dependence and varying temporal dynamics make it difficult to apply the aforementioned works to the theoretical analysis of our time-varying P-values for the robust estimate of the time-varying correlation function.
Our work is also related to locally stationary time series (which models smoothly changing temporal dynamics), see dahlhaus2019towards for a comprehensive review. For theoretical insights, Gaussian approximation proves to be a powerful tool, see zhang2018gaussian for maximum-type Gaussian approximation results of high-dimensional time series, and chernozhukov2017central, chang2024central and wu2024asynchronous for Gaussian approximation of first-order structures of time series in hyperrectangles. In contrast, to analyze numerous P-values from dependent multidimensional inference problems of correlation functions, we derive state-of-the-art Gaussian approximation results for the product sequence of two locally stationary time series in hyperrectangles ($ \{ {\mathbf y} \in \mathbb{R}^p: y_i \le x_i, i \in [p] \} $), where $p$ is allowed to diverge.
Notations. For an integer $N$, denote $[N] = \left\{ 1,2,\dots,N \right\}$. For a vector $\mathbf{v} = (v_1,\cdots,v_p) \in \mathbb{R}^p$, let $|\mathbf{v}| = (\sum_{j=1}^p v^2_j)^{1/2}, |{\mathbf v}|_\infty = \max_j|v_j|$. For another vector $\mathbf{w} = (w_1,\cdots,w_p) \in \mathbb{R}^p$, write ${\mathbf v} \le {\mathbf w}$ if $v_i \le w_i$ for all $i\in[p]$. For a random vector $\mathbf{V}$, $q \geq 1$, let $\|\mathbf V\|_q = (\mathbb{E}|\mathbf V|^q)^{1/q}$. In this paper we consider the kernel function $K(\cdot)$ that is zero outside $(-1, 1)$, and write $K_{b} = K(\cdot/b)$ for some bandwidth $b$. For an index set $B$, let $|B|$ denote its cardinality.
We consider the following time series model:
for dimension $i\in[p]$ and observation $j = [n]$, where $Y_{j,i}$ is the observation, $t_j = j/n$ is the observation time, $(\epsilon_{j,i})_{j=1}^n$ is the locally stationary error process accounting for smooth temporal changes, $\mu_i(\cdot)$ is the deterministic trend function on $[0,1]$ with $d_i$ abrupt change points $0 = a_{i,0} < a_{i,1} < \cdots < a_{i,d_i} < a_{i, d_i+1} =1$. We assume the trend function $\mu_{i}(t)$ is Lipschitz continuous over $[a_{i,m-1}, a_{i, m})$, and the Lipschitz constants are uniformly bounded for $i\in[p],m\in[d_i+1]$.
For each dimension \( i \in [p] \), the local stationarity of \( (\epsilon_{j,i})_{j=1}^n \) indicates that there exists a stationary process \( \{\varepsilon^{\dagger}_{j,i}(u)\}_{j=1}^n \) for each $u\in(0,1)$ such that \(\epsilon_{j,i} = \varepsilon^{\dagger}_{j,i}(j/n)\), and \( \varepsilon^{\dagger}_{j,i}(u) \) satisfies the stochastic Lipschitz condition \(\|\varepsilon^{\dagger}_{j,i}(u) - \varepsilon^{\dagger}_{j,i}(s)\|_q \leq C |u - s|\) for \( u, s \in [0,1] \). This implies that the distribution of \( \epsilon_{j,i} \) changes slowly as \( j \) increases. Local stationarity provides a flexible framework beyond stationarity, and hence has recently attracted significant attention in time series analysis.
We study the case $d_i > 0$. The cases where $d_i = 0$ can be analyzed similarly.
For dimensions $i,l \in [p]$, we define the following time-varying covariance, marginal variance and correlation function as
Here, $\gamma_{i,l}(t)$ is a function defined on $[0,1]$, independent of $n$, and can be approximated by the covariance between the discrete-time processes $\epsilon_{\lfloor nt \rfloor, i}$ and $\epsilon_{\lfloor nt \rfloor, l}$ as $n \to \infty$. It will become constant if the error process is stationary. We assume $\gamma_{i,l}(t)$ is twice continuously differentiable, which is proposed by many works, see zhao2015inference, CHEN2025network for examples.
The time-varying correlation network is defined as ${\cal G}(t) = \left\{ {\cal V},{\cal E}(t) \right\}$ with vertex set ${\cal V} = [p]$, and edge set ${\cal E}(t) = \left\{ (i,l)\in {\cal V}\times{\cal V}: \gamma_{i,l}(t) \ne 0, i\ne l \right\}$. Since $\gamma_{i,l}(t) = \gamma_{l,i}(t)$, learning the correlation network at time $t\in[0,1]$ is equivalent to testing the following $p(p-1)/2$ hypotheses: for $(i,l)\in{\cal H}$,
where ${\cal H} = \left\{ (i,l):i,l\in[p],i > l \right\}$. An edge in the network is established between node $i$ and node $l$ at time $t$ if and only if $H_{i,l,t}^0$ is rejected.
Let ${\bf Y}$ denote the collection of all observations, ${\cal R}({\bf Y},t) = \{(i,j)\in{\cal H}: H_{i,l,t}^0 ~\mbox{is rejected}\}\subset {\cal H}$ be the constructed edge set at time $t$, and ${\cal H}_0(t)= \left\{ (i,l)\subset{\cal H}: \rho_{i,l}(t)=0 \right\}$ denote the true edge set at time $t$. The false discovery proportion (FDP) is defined as
and the FDR is defined as
In our context, the FDP is defined as the ratio of false connections in our correlation network to the total number of connections identified, and the FDR represents the expected value of this ratio. Our objective is to control FDR under the specified level $\alpha$ for all $t\in(0,1)$.
To estimate the correlation function $\rho_{i,l}(t)$ in the presence of abrupt change points in means, the difference-based approach provides an effective solution, as illustrated in bai2025time. Prefix a sequence $h\rightarrow\infty, h = o(n)$. For $j\in[n],i,l\in[p]$, define $\beta_{i,l}(t_{j}) := \mathbb{E}\left\{ (\epsilon_{j,i} - \epsilon_{j-h,i})(\epsilon_{j,l} - \epsilon_{j-h,l}) \right\}$. If $h$ is large enough, then $\mathbb{E}(\epsilon_{j,i}\epsilon_{j-h,l})\approx 0$, and simple calculations will give $\beta_{i,l}(t_{j}) \approx 2\gamma_{i,l}(t_j)$. Additionally, if the trend is relatively smooth, we have $y_{j,i} := Y_{j,i}-Y_{j-h,i}\approx \epsilon_{j,i} - \epsilon_{j-h,i}$. Thus we can estimate $\gamma_{i,l}(t)$ by the local linear estimator $\tilde \gamma_{i,l}(t) = \hat \beta_{i,l}(t)/2$ satisfying
where $b_{i,l}$ is the bandwidth parameter. In the notation, we suppress the dependence of all these on $h$, and the choice of $h$ is deferred to Section (ref) of the Appendix in for discussion. Consequently the cross-correlation function $\rho_{i,l}(t)$ can be estimated by
We establish a general framework to obtain time-varying P-values which involves two key steps: first, expand the maximum deviation of correlation estimation by a weighted sum of locally stationary processes, applying Gaussian approximation theory and multiplier bootstrap to obtain bootstrap samples that can approximate the distribution of maximum deviation; second, based on the bootstrap samples, obtain a set of time-varying P-values for further FDR control.
Define $e_{j,i,l} = (\epsilon_{j+h,i} - \epsilon_{j,i})(\epsilon_{j+h,l} - \epsilon_{j,l}) -\beta_{i,l}(t_{j})$. For $i,l\in[p]$, from (ref) and (ref) in the Appendix we could approximate the maximum deviation of the cross-correlation estimate via $\vartheta_{i,l}(t)$, the moving weighted average of innovations, i.e.,
and for each $i,l\in[p]$, $(\Xi_{j,i,l})_{j=1}^n$ is a locally stationary innovation process. Let $c_{i,l} = (b/b_{i,l})^{1/2}$, where $b=\max_{i,l\in[p]} b_{i,l}$. Define the long run variance of $\vartheta_{i,l}(t)$ and the standardized version of $(\Xi_{j,i,l})_{j=1}^n$ as:
Using the decomposition of (ref) and (ref), we could use the following approximation intuitively:
For the estimation of $\Xi_{j,i,l}$, $\Gamma_{i,t}^2(t)$ $\bar{ \Xi}_{j,s,i,l}$, we follow the same method as in bai2025time to estimate them with $\tilde{\Xi}_{j,i,l}$, $\tilde{\Gamma}_{i,t}^2(t)$ and $\hat{\bar{ \Xi}}_{j,s,i,l}$; the details are postponed to Section (ref) in the Supplementary.
The multiplier bootstrap in step 4-7 of Algorithm (ref) is justified in Theorem (ref). In Theorem (ref) we show that our P-values in (ref) are simultaneously and asymptotically valid for all $t\in[\max_{i,l}b_{i,l},1-\max_{i,l}b_{i,l}]$ and exhibit high power in rejecting false null hypotheses.
Let $P_{(1)}(t)\le \dots\le P_{(p(p-1)/2)}(t)$ denote the order statistics of the P-values at time $t$. Define a sequence of increasing threshold $0\le\Delta(1)\le\dots\le \Delta\{p(p-1)/2\} \le 1$. At time $t$ we reject the hypotheses with P-values smaller than $\Delta\left\{ R({\bf Y},t) \right\}$, where $R({\bf Y},t) = \max\left\{ r:P_{(r)}(t) \le \Delta(r) \right\}$ is the number of rejected hypotheses. To control FDR level under $\alpha$ the Benjamini--Hochberg (B-H) procedure takes
and Benjamini--Yekutieli (B-Y) procedure takes
where $\gamma$ is Euler's constant.
The B-H procedure effectively controls the FDR when P-values are independent or satisfy the positive regression dependence on a subset (PRDS) benjamini2001control, see also our assumption (ref). For arbitrarily dependent P-values, the B-Y procedure is valid benjamini2001control, although it is known to be conservative.
When implementing Algorithm (ref), we need to select the following parameters: the lag $h$, the bandwidths $b_{i,l}$ used for estimating $\rho_{i,l}(t)$, the parameters $m_{i,l}$ and $\eta$ for estimating $\Gamma_{i,l}(t)$, as well as the window length $w$ in Step 3 of Algorithm (ref). We outline our methodology in Section (ref) of the Appendix. The simulation results in Sections (ref) and (ref) of the Appendix demonstrate the effectiveness and robustness of our approach.
We assume for each dimension $i\in[p]$, the error process $(\varepsilon_{j,i})_{j\in[n]}$ admits $$ \varepsilon_{j,i} = G_i(t_j,{\cal F}_j),$$ where $t_j=j/n$, $G_i(\cdot,\cdot)$ are measurable functions of $[0,1]\times \mathbb{R}^{\mathbb{Z}}\rightarrow\mathbb{R}$, ${\cal F}_j= (\cdots,\xi_{j-1},\xi_j)$, and $(\xi_j)_{j\in\mathbb{Z}}$ are i.i.d. random elements. This modeling of $\varepsilon_{j,i}$ is very flexible and includes the time-varying autoregressive model and time-varying ARCH and GARCH model (see Dannazhang2021). For any process $G(t,{\cal F}_j)$, we say it is ${\cal L}^q$ stochastic Lipschitz continuous with constant $C$ (denoted as $G(\cdot,\cdot)\in\mathrm{Lip}_q$), if for $t_1,t_2\in[0,1]$, $\left\| G(t_1,{\cal F}_0)-G(t_2,{\cal F}_0) \right\|_q\le C|t_1-t_2|$. For $q \ge 2$, such ${\cal L}^q$ stochastic Lipschitz continuous process is called a locally stationary process. For $j \ge 0$, let ${\cal F}_j^* = ({\cal F}_{-1},\xi_0',\xi_1,\dots,\xi_j)$, where $\xi_0'$ is an i.i.d. copy of $\xi_0$. The physical dependence measure for $G(\cdot,\cdot)\in\mathrm{Lip}_q$ is defined by $$\delta_q(G,l) = \sup_{t\in[0,1]}\left\| G(t,{\cal F}_l)-G(t,{\cal F}_l^*) \right\|_q$$ if $l\ge 0$ and $\delta_q(G,l) = 0$ if $l < 0$. Physical dependence measure quantifies the influence of the input $\varepsilon_0$ on the output $G(t,{\cal F}_l)$, which is an effective measure of the strength of dependence.
Based on the dependence measure above, we propose the following assumptions about moments, dimensionality, dependence strength, regularity of the error process, structural properties of trend functions with jumps, kernels and bandwidths.
Assumptions (ref)(\romannumeral 1) and (ref)(\romannumeral 3) are standard in the kernel-based nonparametric analysis of locally stationary time series, see Assumption 1 of zhao2015inference for examples. Assumption (ref)(\romannumeral 2) indicates sub-Weibull tails (weaker than sub-exponential conditions). Assumption (ref)(\romannumeral 4) and (ref)(\romannumeral 5) guarantees $\gamma_{i,l}(t), \Gamma_{i,l}(t)$ are smoothly changing with respect to $t$. Assumption (ref)(\romannumeral 6) accommodates scenarios with an increasing number of change points. Assumptions (ref)(\romannumeral 7) and (ref)(\romannumeral 8) are standard regularity conditions related to the kernel and bandwidth.
Recall that $w$ is the block size in Step 3 of Algorithm (ref), $m_{i,l}$ and $\eta$ are the smoothing parameters for the estimator of $\Gamma_{i,l}(\cdot)$, see (ref). Let $m = \max_{i,l\in[p]} m_{i,l}$, $\vartheta_n = \frac{\log^2 n}{w} + \frac{w}{nb} + (\frac{w}{nb})^{1/2}(np^2)^{4/q}$, $g_n = p^{2/q}\big\{(\frac{m}{n \eta^{2}})^{1/2}+m^{-1}+\eta +(\frac{m}{nb})^{1/2}(\frac{mb}{n})^{-1/(2q)}+ (nb)^{-1/2}b^{-1/q}\big\} + w^{3/2}/n + n^{\phi}h/\sqrt{w}$.
Assumption (ref) can be satisfied for sufficiently large $q$ if $\omega\asymp n^{2/5},\eta\asymp n^{-1/7}, b\asymp n^{-1/5}, m\asymp n^{2/7}, \phi < 1/5$ and $p\asymp n^{\iota}$, where $\iota > 0$ is a fixed positive number. It allows $p$ to increase polynomially with respect to sample size $n$.
Intuitively, the PRDS condition means that rejecting hypothesis $(i, l)$ at time $t$ increases the probability that other hypotheses are rejected at time $t$. In other words, this condition reflects a form of positive dependence among the P-values.
Recall the definition of $\bar { \Xi}_{j,s, i,l}$ in (ref). For each $i,l\in[p], i >l$, define an $(n-2\lceil n b \rceil + 1)$-dimensional vector $\bar{\mathbf \Xi}_{j,i,l} = (\bar{ \Xi}_{j,\lceil n b \rceil,i,l}, \bar{\Xi}_{j + 1,\lceil n b \rceil + 1, i, l}, \cdots, \bar{\Xi}_{n - 2\lceil n b \rceil +j, n - \lceil n b \rceil,i,l})^{{ \mathrm{\scriptscriptstyle T} }}$, and $\bar{\mathbf \Xi}_{j}^{{\cal G}} = (\bar{\mathbf \Xi}_{j,2,1},\bar{\mathbf \Xi}_{j,3,1},\bar{\mathbf \Xi}_{j,3,2},\dots,\bar{\mathbf \Xi}_{j,p,p-1})^{{ \mathrm{\scriptscriptstyle T} }}$ as an $(n-2\lceil n b \rceil + 1)p(p-1)/2$-dimensional vector. To simplify notation, for any $np(p-1)/2$-dimensional vector ${\bf Z} = ({\bf Z}_{2,1},{\bf Z}_{3,1},{\bf Z}_{3,2},\dots,{\bf Z}_{p,p-1})^{ \mathrm{\scriptscriptstyle T} }$ with ${\bf Z}_{i,l} = (Z_{1,i,l},\dots, Z_{n,i,l})^{ \mathrm{\scriptscriptstyle T} }$, define
The Gaussian approximation in bai2025time is of the maximum-type, which is a special case of (ref) with ${\mathbf x} = (a,\dots,a),a\in\mathbb{R}$, and it suffices for FWER control. However, FDR control requires a Gaussian approximation on hyper-rectangles (derived from chang2024central and wu2024asynchronous, with necessary modifications to adapt to our setting) to ensure valid inference and P-value computation in each dimension.
Theorem (ref) and (ref) imply that, the bootstrap samples in the steps 4--7 of Algorithm (ref) can effectively approximate the maximum deviation of the standardized cross-correlation estimates. Based on this, in Theorem (ref), we show that our P-values $P_{i,l}(t)$'s (defined in (ref)) are asymptotically valid for any time $t\in{\cal T} = [b,1-b]$, and have sufficient power to reject the null hypothesis whenever $\rho_{i,l}(t)\ne 0$.
Theorem (ref) and provides theoretical guarantees for the B-H and B-Y procedure. Recall the definition of FDR in (ref).
In this section, we present the results of our simulations. We apply Algorithm (ref) in conjunction with the FDR control procedures in Section (ref). The performance of our approach is evaluated using the following two metrics: (\romannumeral 1) FDP control: The FDP should be controlled below the prescribed level. (\romannumeral 2) Low FNP: The FNP of our method should decrease rapidly as the sample size increases, demonstrating the method’s effectiveness in detecting true signals.
Besides, the results should not exhibit abrupt changes, as the difference-based approach circumvents the need for pre-estimating change points. Finally, the results should accurately capture and reflect the time-varying nature of the underlying process.
Recall the model ${\bf Y}_j=\boldsymbol \mu(t_j)+{\bf G}(t_j,{\cal F}_j), j\in[n]$. In the following, we present details of our simulation setting.
Case 1 ($p = 6$): Let $\left\{ \eta_{j,i},j\in[n],i\in[6] \right\}$ be the collection of i.i.d.\ standard Gaussian random variables and define $\boldsymbol{\eta}_j = (\eta_{j,1}, \dots, \eta_{j,6})^\top$. Define the function $f_0(t) = 0.25 - 0.1 (t - 1/2)^2$. Let ${\bf I}_p$ denote the $p$-dimensional identity matrix, ${\bf J}_p$ denote $p$-dimensional matrix with each entry being $1$, and ${\bf A}(t) = f(t){\bf I}_6$. For the error process, we set $\mathbf{G}(t, \mathcal{F}_j) = \mathbf{A}(t) \mathbf{G}(t, \mathcal{F}_{j-1}) + (4{\bf I}_6/5 + {\bf I}_2\otimes{\bf J}_3/5) \boldsymbol{\eta}_j$. The mean functions are defined as
where for $j = 1,4$, $a_{j,1} = 0.25$ and $a_{j,2} = 0.55$; for $j = 2,5$, $a_{j,1} = 0.4$ and $a_{j,2} = 0.7$; and for $j = 3,6$, $a_{j,1} = 0.55$ and $a_{j,2} = 0.85$. There are $p(p-1)/2 = 15$ total hypotheses with $|\mathcal{H}_0| = 9$ true nulls, resulting in a target FDR of $0.06$.
Case 2 ($p = 9$): Similar to Case 1, let $\left\{ \eta_{j,i},j\in[n],i\in[9] \right\}$ be a collection of i.i.d.\ standard Gaussian random variables and define $\boldsymbol{\eta}_j = (\eta_{j,1}, \dots, \eta_{j,9})^\top$. The error process is $\mathbf{G}(t, \mathcal{F}_j) = \mathbf{A}(t) \mathbf{G}(t, \mathcal{F}_{j-1}) + (4{\bf I}_9/5 + {\bf I}_3\otimes{\bf J}_3/5) \boldsymbol{\eta}_j$, where ${\bf A}(t) = f(t){\bf I}_9$. For the mean function, let $\mu_i(t)$, $i \in [6]$, be the same as in Case 1, and set $\mu_7(t) = \mu_1(t)$, $\mu_8(t) = \mu_2(t)$, and $\mu_9(t) = \mu_3(t)$. There are $p(p-1)/2 = 36$ total hypotheses with $|\mathcal{H}_0| = 27$ true nulls, yielding a target FDR of $0.075$.
We consider sample sizes $n \in \{450, 600\}$ with lag choices $h = \lceil 2 \log n \rceil$. We apply Algorithm (ref) followed by the B-H procedure at level $\alpha = 0.1$ for Cases 1 and 2. Additionally, we apply the B-Y procedure at level $\alpha = 0.1$ for Case 1 for comparison. Each experimental setting was repeated 100 times, and the reported results are averages over these repetitions.
Table (ref) presents the maximum FDP and FNP over $t \in (1/3, 2/3)$ (peak proportions), as well as the mean FDP and FNP over the same interval (average proportions). Figure (ref) are the time-varying FDP and FNP trajectories over $t \in (1/3, 2/3)$. For simplicity, only the trajectories of Case 1 with the B-H procedure are presented; trajectories for other cases and methods exhibit similar patterns.
Based on the above results, we made the following observations:
(\romannumeral 1) FDP and FNP levels: As shown in Table (ref), our method effectively controls the maximum FDP close to the preassigned levels while maintaining a relatively low FNP level.
(\romannumeral 2) Change point issue: As shown in Figure (ref), the smooth trajectories of the FDP and FNP indicate that our method successfully handles jump points in observations.
(\romannumeral 3) Time-varying nature: The FNP and FDP levels in Figure (ref) vary over time, indicating that our method can effectively capture the time-varying nature of the data.
(\romannumeral 4) B-H and B-Y methods: As shown in Table (ref), the B-H method yields lower FNPs and FDPs closer to the target level, while the B-Y method achieves lower FDPs but with a higher FNPs. In both methods, FDPs are well controlled and FNPs decrease rapidly as sample size increases.
We also compare our method to the moving-window approach described in Section 3.9 of masuda2025introduction. This method selects a window length $w$ and, at time $t = j/n$, computes sample correlations using data from the time interval $[j-w, j+w]$. Two nodes are connected at time $t$ if their sample correlation exceeds a chosen threshold. As pointed out in masuda2025introduction, there is no consensus on an optimal method for selecting the threshold values. Therefore, in Table (ref) we present the results of moving-window approach across several different threshold levels, using Case 1 as a representative example. The results show that, depending on the threshold, the method either produces an excessively high FDP or is overly conservative, highlighting the challenges of threshold selection in this approach.
Section (ref) of the Appendix deals with sensitivity analysis, showing that our method is robust to slight perturbations in the bandwidth and lag parameters. Section (ref) of the Appendix compares the performance of our method with the method in bai2025time, which focuses on the FWER control.
Electroencephalography (EEG) measures the electrical activity of the brain using sensors placed on the scalp and facilitates many neuroscience studies. In this section, we analyze an EEG dataset collected to explore the relationship between alcoholism and brain activity. Details regarding data collection procedures can be found in zhang1997electrophysiological. While the original analysis focuses on amplitudes of event-related potentials (ERPs), we adopt an alternative perspective by constructing time-varying correlation networks of the ERPs.
The data in our analysis are available at \url{https://archive.ics.uci.edu/dataset/121/eeg+database}, and we focus on the SMNI_CMI_TRAIN dataset. It comprises 10 subjects from the alcoholic group and 10 subjects from the control group. Each subject performed 10 image recognition tasks when a single object is shown, with each task lasting one second, resulting in 10 segments per subject. Following the methodology in zhang1997electrophysiological, the ERPs for each subject were averaged across these segments to obtain mean ERPs. ERPs were recorded using $p = 64$ electrodes, and each ERP segment consists of $n = 256$ time points, corresponding to a sampling rate of 256 Hz over one second. For each subject, we construct the time-varying correlation network $\mathcal{G}(t) = \{\mathcal{V}, \mathcal{E}(t)\}$ based on mean ERPs using Algorithm (ref) in conjunction with the B-H procedure. For comparison, we also analyze the data with the FWER control method described in bai2025time.
The 64 scalp electrodes can be grouped into five anatomical regions: frontal, central, parietal, occipital, and temporal. Let $\mathcal{V}_1, \dots, \mathcal{V}_5$ denote the sets of nodes corresponding to these five anatomical regions, and ${\cal V}$ denote the set of all nodes. We define the connection proportion for region $i$ at time $t$ as
Figure (ref) displays the time-varying connection proportions for each region (averaged across 10 subjects in each group). We have the following observatios:
(\romannumeral 1) The connection proportions for each region exhibit clear temporal variations, reflecting changes in brain activity in response to the image stimuli presented to the subjects.
(\romannumeral 2) The control subjects demonstrate significantly higher level of correlations, in the temporal, central, and frontal regions, and a slightly higher level of correlation in the occipital region. Similarly, zhang1997electrophysiological reported stronger brain activity in the temporal and occipital regions of control subjects. While their analysis focused on amplitudes, our results highlight differences in correlation as well.
(\romannumeral 3) The results from the FWER control show that the proportions of connections in the temporal and occipital region are very low for both groups. As a result, we cannot conclude that correlations are higher in these regions based solely on the FWER control results, despite strong indications from the FDR control method and previous findings reported by zhang1997electrophysiological. This suggests that the FWER approach may be excessively stringent, potentially masking some meaningful results.
The interconnectedness among economic data plays a crucial role in quantifying risks within financial portfolios, as well as systemic and macroeconomic risks. For an introduction, see acemoglu2012network, puliga2014credit, and guo2018development, among others. In this section, we construct time-varying correlation networks of Daily WRDS World Indices, which are market-capitalization-weighted indices with dividends at daily frequency (for details, see \url{https://wrds-www.wharton.upenn.edu/pages/get-data/world-indices-wrds/}).
Specifically, we analyze the weekly averages of the absolute values of the Daily Country Return with Dividends (PORTRET) from June 1, 2006, to June 1, 2022, which can serve as measures of volatility and risk. We denote the weekly average of the absolute value of PORTRET at time \( j \) for region \( i \) by \( Y_{ji} \). There are $n=819$ time points for $p = 6$ regions encompassing France, Germany, the UK, China, Japan, and Hong Kong. Based on these observations, we apply Algorithm (ref) together with the B-H and B-Y procedures, using a significance level of \( \alpha = 0.05 \). The results are presented in Figure (ref). For brevity, the figure displays only the time-varying networks corresponding to the 20th week of each year. From the figure, we highlight the following observations:
(\romannumeral 1) The network construction results accurately reflect geopolitical realities. For example, using Algorithm (ref) with the B-H method, we observe persistent connections over time between France and Germany, Hong Kong and China, China and Japan, and the UK and Germany, consistent with their strong economic relationships.
(\romannumeral 2) The results reflect the temporal variability of correlations. For instance, the connectivity between Hong Kong and other regions exhibits two notable declines: one between 2012 and 2013, and another in 2018. These periods correspond to stock market crashes in Hong Kong, characterized by sharp increases in volatility that eroded confidence in cross-regional correlations. Additionally, prior to 2013, no significant connection existed between China and France; however, in 2014, following the signing of a comprehensive strategic partnership agreement aimed at increasing bilateral trade to 100 billion euros by 2015, our network detects a corresponding emergence of ties.
(\romannumeral 3) When the cost associated with detecting spurious connections is considered too high, researchers may prefer a sparser and more conservative network. In such cases, the B-Y method provides a favorable alternative. We also present the results obtained via the B-Y method in Figure (ref), where the network is sparser.
To demonstrate the effectiveness of our method in capturing correlations, we additionally construct semi-synthetic data. We introduce a synthesized country whose the weekly average of the absolute value of PORTRET at time $t$ is defined as follows. Let $\widetilde{Y}_j,j\in[819]$ be obtained by randomly sampling without replacement from $\left\{ Y_{ij},i\in[6],j\in[819] \right\}$, and $Y_{j,1}$ represents the data of China. We costruct the synthesized country data as $\left\{ Y_{j,7},j\in[819] \right\}$, where
The correlation signal emerges at \( t = 1/3 \) and is effectively detected at \( t = 0.406 \) by our Algorithm (ref) with B-H method at significance level $\alpha = 0.05$. This result on semi-synthesized data demonstrates that our method is capable of detecting signals accurately and efficiently.
In this paper, we propose a novel framework for learning time-varying correlation networks that adapt to complexity of data, including high-dimensionality, non-Gaussianity, non-stationarity and abrupt change points in the mean, while effectively controlling FDR. We first employ a difference-based approach to estimate time-varying correlation functions, mitigating the influence of abrupt change points. Then we develop an advanced Gaussian approximation theory and generate novel time-varying P-values via bootstrap to test for the presence of time-varying correlations. Furthermore, we combine these time-varying P-values with the B-H and B-Y procedures, providing theoretical guarantees for FDR control allowing for cross-sectional dependence. The effectiveness of our methodology is validated through extensive simulations. In addition, we demonstrate the practical utility of our method by analyzing both EEG and financial datasets, yielding meaningful and interpretable results.
This work has several limitations, but it also paves the way for some exciting future research. First, our study primarily focuses on correlations; however, partial correlations also play a crucial role in characterizing the interactions between variables and deserve consideration in future work. Additionally, we have limited our investigation to the well-known B-H and B-Y procedures for combining time-varying P-values. Exploring alternative methods for FDR control could yield further insights. Lastly, the PRDS condition adopted in this work may be restrictive for certain applications. Therefore, it would be valuable to examine other dependence structures in FDR control frameworks in future work.
We gratefully thank Prof. Huaqing Jin for providing background information and some relevant references on brain EEG data. Weichi Wu is the corresponding author. Bufan Li and Weichi Wu are supported by the NSFC 12271287. Lujia Bai is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation, TRR 391: Spatio-Temporal Statistics for the Transition of Energy and Transport under Project 520388526 and Research unit 5381 Mathematical Statistics in the Information Age, project number 460867398).