Online Correlation Change Detection for Large-Dimensional Data with An Application to Forecasting of El Niño Events

Consider a sequence of observations $\{\bm{x}_{t},t=1,2,\ldots\}$, where $\bm{x}_{t}\in{\mathbb{R}}^{p}$. In general, we do not impose any restriction on $p$, except in the knockoff enhancement discussed in Section 3.4. We are interested in detecting the change point in the underlying correlation structure. We assume the following data-generating model:

Here, $\nu$ is the change-point at which the correlation structure of $\bm{x}_{t}$ changes. The matrices $\bm{R}_{0}=[\rho_{0}(i,j)]_{i,j=1,\ldots,p}$ and $\bm{R}_{1}=[\rho_{1}(i,j)]_{i,j=1,\ldots,p}$ are the pre- and post-change correlation matrices, respectively. We assume that both $\bm{R}_{0}$ and $\bm{R}_{1}$ are unknown, and the change point $\nu$ is deterministic and unknown. The difference between $\bm{R}_{0}$ and $\bm{R}_{1}$ indicates the magnitude and pattern of the change. Specifically, the correlation change is considered sparse if the number of differing entries between $\bm{R}_{0}$ and $\bm{R}_{1}$ is relatively small compared to the total number of entries, otherwise the change is considered dense.

Our goal is to detect the change as quickly as possible, subject to false alarm constraints. This is usually achieved by designing a stopping time (Poor and Hadjiliadis, 2008), which is a random variable $T$ relating to the data sequence $\{\bm{x}_{t},t=1,2,\ldots\}$ such that for each $t$, the event $\{T=t\}\in\sigma(\bm{x}_{1},\ldots,\bm{x}_{t})$, where $\sigma(\bm{x}_{1},\ldots,\bm{x}_{t})$ denotes the sigma-algebra generated by $\{\bm{x}_{1},\ldots,\bm{x}_{t}\}$. Equivalently, the event $\{T=t\}$ is a function of only $\bm{x}_{1},\ldots,\bm{x}_{t}$.

A fundamental objective in sequential change point detection is to optimize the trade-off between false alarm rate and average detection delay. Controlling a false alarm rate is commonly achieved by setting an appropriate threshold on a test statistic. On the other hand, the threshold also affects the average detection delay. A larger threshold incurs fewer false alarms but leads to a larger detection delay, and vice versa.

We introduce next the two commonly used performance metrics in sequential detection, the average run length (ARL) and the worst-case average detection delay (WADD). ARL is used to characterize false alarms and it is defined, for a given stopping time $T$, as:

its reciprocal is the commonly used false-alarm rate ${\mathsf{FAR}}(T)=1/{{\mathsf{ARL}}(T)}=1/{\mathbb{E}_{\infty}[T]}$. ARL can be interpreted as the expected time duration between two consecutive false alarms. We typically focus on the test procedures that satisfy a constraint $\gamma$ on the ARL, i.e., consider the set of tests:

notice that ARL (or $\gamma$) is pre-specified and can go to infinity.

Finding a uniformly powerful test that minimizes the delay over all possible values of the change point $\nu$, subject to an ARL constraint, is generally intractable. Usually, it is more tractable to pose the problem in the so-called minimax setting. We adopt the worst-case detection delay (WADD), defined as the supremum of the average detection delay conditioned on the worst possible realizations (Lorden, 1971). More specifically,

where $\operatorname*{ess\,sup}$ is essential supremum, and $(x)^{+}=\max\{x,0\}$. Under such a definition of WADD, the formulation of interest is:

Assume a data sequence $\bm{x}_{1},\ldots,\bm{x}_{t},\ldots$ that is i.i.d., where $\bm{x}_{t}:=(x_{t1},\ldots,x_{tp})^{\top}$ with bounded fourth-moment entries and correlation matrix $\bm{R}$. The matrix $\bm{R}$ can be estimated using the sample correlation matrix defined as follows. There are typically two scenarios, depending on the prior knowledge of the distributional parameters. For illustration purposes, we first consider the simplest case where the mean and variance of each entry are known. Without loss of generality, we assume $\mathbb{E}(x_{1i})=0$ and $\sigma^{2}=\mathbb{V}\mathrm{ar}(x_{1i})=1$, $\forall i=1,\ldots,p$. The sample correlation matrix for samples within a time window $[s,t]$ is then given by $\bm{\widetilde{R}}_{s:t}=\frac{1}{t-s+1}\sum_{k=s}^{t}\bm{x}_{k}\bm{x}_{k}^{\top}$.

Second, we consider the general case where the population mean and variance are unknown. Let $\bm{\bar{x}}_{s:t}=\frac{1}{t-s+1}\sum_{k=s}^{t}\bm{x}_{k}$ be the sample mean from time $s$ to $t$, and the corresponding sample covariance matrix within $[s,t]$ is

Let $\bm{D}_{s:t}=diag(\bm{\widehat{\Sigma}}_{s:t})$ be the diagonal matrix of $\bm{\widehat{\Sigma}}_{s:t}$ and $\bm{y}_{k}=\bm{D_{s:t}}^{-1/2}(\bm{x}_{k}-\bm{\bar{x}}_{s:t})$ be the standardized vector of $\bm{x}_{k}$, for $k=s,\ldots,t$, then the corresponding sample correlation matrix for samples within $[s,t]$ is

Assume there is also a set of reference data (historical data) $\{\bm{x}_{-H},\ldots,\bm{x}_{0}\}$ of size $H+1$, which are known to be generated from the pre-change distribution with an unknown correlation matrix $\bm{R}_{0}$. According to equation (5), the sample correlation matrix of all reference data can be similarly calculated as:

where $\bm{y}_{k}=\bm{D}_{-H:0}^{-1/2}(\bm{x}_{k}-\bm{\bar{x}}_{-H:0})$, for $k=-H,\ldots,0$, and we further denote $\bm{\widehat{R}}_{0}=[\hat{\rho}_{0}(i,j)]_{i,j=1,\ldots,p}$. The size of historical data, $H+1$, is assumed to be sufficiently large. When appropriate, we may use only a subset of this historical data for calculating the sample correlation matrix.

Given the online sequence $\{\bm{x}_{t},t=1,2,\ldots\}$, we calculate the detection statistics as follows. At each time $t$, we investigate all values of $t^{\prime}<t$ as potential change-points. For each $t^{\prime}$, we calculate the sample correlation matrix $\bm{\widehat{R}}_{t^{\prime}:t}=[\hat{\rho}_{t^{\prime}:t}(i,j)]_{i,j=1,\ldots,p}$ according to (5) using potential post-change samples $\{\bm{x}_{t^{\prime}},\ldots,\bm{x}_{t}\}$. Then, we calculate the squared difference between $\bm{\widehat{R}}_{t^{\prime}:t}$ and $\bm{\widehat{R}}_{0}$, i.e.,

where ${vecho}(\cdot)$ indicates the half-vectorized vector by vectorizing only the lower triangular part (excluding the diagonal) of the symmetric matrix. If a change point exists at $t^{\prime}$, then $\bm{v}_{t^{\prime},t}$ is obviously an estimator of $vecho\big{\{}[\left(\rho_{0}(i,j)-\rho_{1}(i,j)\right)^{2}]_{i,j=1,\ldots,p}\big{\}}$, which measures the entry-wise differences between the pre- and post-change population correlation matrices.

More specifically, for any $(i,j)$-th entry of the correlation matrix, we present the expectation of $\bm{v}_{1,t}(i,j)$ when the data sequence $\{\bm{x}_{t},t=1,2,\ldots\}$ is sampled entirely from the post-change regime, to illustrate its effectiveness for detection. For notational convenience, we write $\rho_{0}=\rho_{0}(i,j)$ and $\rho_{1}=\rho_{1}(i,j)$ here.

We note that in (8), the first term $(\rho_{0}-\rho_{1})^{2}$ dominates the expectation when $H$ and $t$ are sufficiently large. Therefore, for any fixed $t$, larger values of $\bm{v}_{t^{\prime},t}(i,j)$ indicate a significant difference between $\rho_{0}(i,j)$ and $\rho_{1}(i,j)$. In contrast, under the pre-change regime, all elements in the $\bm{v}_{t^{\prime},t}$ are considerably small.

To detect dense changes in correlation structure, we use the $\ell_{1}$ norm of $\bm{v}_{t^{\prime},t}$ to construct detection statistics. Since the true change-point location $\nu$ is unknown, we take the maximum for all potential change-point $t^{\prime}<t$, and the resulting detection statistic is defined as

the scaling $\frac{(t-t^{\prime})H}{H+t-t^{\prime}}$ is incorporated to balance the variance of $\bm{v}_{t^{\prime},t}$, as the variance of $\bm{v}_{t^{\prime},t}$ is $O(\frac{1}{t-t^{\prime}}+\frac{1}{H})$. When $H$ is relatively large, we may simply use $S_{t}^{(\text{sum})}=\max_{1\leq t^{\prime}\leq t-1}(t-t^{\prime})\|\bm{v}_{t^{\prime},t}\|_{1}$.

The sum-type detection statistics above may not be efficient in detecting sparse changes, where only a relatively small number of elements change in the correlation structure. Under the sparse setting, the conventional test statistic is the maximum difference between the pre- and post-change correlation matrices. We thus define the max-type test statistic as

For window-limited-variant test statistics, the underlying change point $t^{\prime}$ can vary from $t-w$ to $t-1$, the calculations of $\bm{\widehat{R}}_{t^{\prime}:t}$ and statistics can be unstable due to lack of samples, especially when $w$ is small compared to the dimension of data $p$ and $t^{\prime}$ is close to $t-1$. To solve this issue, we come up with two modified approaches combining with the SMOTE and Knockoff techniques, respectively. See Figure 1 for an illustration.

Selecting a proper threshold $b$ for stopping time is crucial in the sequential change point detection procedure. In general, there are two types of methods for determining the appropriate threshold. One way is to infer an exact threshold based on the distribution of test statistics under the pre-change regime. The other way is to determine the threshold based on the empirical distribution of numerous simulated test statistics, which is more common in practice as it requires little knowledge about the distribution of test statistics.

However, in real-data applications where the number of pre-change samples is limited, or when generating samples is time-consuming, the available data size is small. We propose a new heuristic method for selecting the threshold using a sign-flip permutation method, which works well even for small sample sizes. The method is detailed in Algorithm 2. It is worthwhile mentioning that we assume that the prechange correlation matrix $\bm{R}_{0}$ is known or preestimated using available historical data in Algorithm 2 for simplicity. In real applications, $\bm{{R}}_{0}$ may be known in advance based on domain knowledge or estimateable from reference data, as shown in our real data analyses in Section 6.

To illustrate why the threshold selected via Algorithm 2 is valid, we present the following lemma.

Lemma 3.7 demonstrates that for a given vector $\bm{x}_{t}$, the data shuffled through various sign-flip steps using different $\bm{\mathcal{R}}$ exhibit weak dependence. Moreover, if the shuffled data follows a Gaussian distribution, it becomes independent. This property ensures that the test statistics computed from $q$ sign-flip trials are nearly independent, allowing for appropriate threshold selection. A numerical result supporting this assertion is presented in Remark 3.8.

For simplicity, we let the size of the reference data $H$ equal the chosen window size $w$. Note that we can regard the weight $\frac{(t-t^{\prime})H}{H+t-t^{\prime}}$ before the test statistic $S_{t}$ as a constant for a given $t^{\prime}$, and we may unify both sum-type statistics for non-sparse settings and max-type statistics for sparse settings since they are, in essence, the combination of element-wise entries. We first characterize the temporal correlation of the statistics as follows, which will be utilized later for the ARL approximation.

Based on this temporal correlation, we have the following Theorem characterizing the approximate ARL of the proposed detection procedures. The main idea is to use a linear approximation for the correlation between detection statistics $S_{t}$ and $S_{t+s}$. Then, the behavior of the detection procedure can be related to a random field. By leveraging the localization theorem (Siegmund et al., 2010), we can obtain an asymptotic approximation for ARL when the threshold $b$ is large enough.

The bounded fourth-moment assumption is needed to ensure that the estimated correlation coefficients are reliable estimators, i.e., $\hat{\rho}(i,j)$ converge to $\rho(i,j)$ in probability. Additionally, the Assumption C.1, C.2, C.3 and C.4 are required to ensure the applicability of the localization theorem. In practice, all assumptions can be easily satisfied by data from some commonly used distributions, such as Gaussian distributions and exponential families.

The main contribution of Theorem 4.2 is to provide a theoretical method to set the threshold that can avoid the Monte Carlo simulation, which could be time-consuming, especially when the desired ARL is large. From Theorem 4.2, we can numerically compute the threshold value $b$ by setting the right-hand side of Equation (15) to the desired ARL value. Table 3 shows the high accuracy of this approximation result by comparing the threshold obtained from Equation (15) with that obtained from a simulation study. For a sequence of i.i.d standard normal samples, we conduct 1000 sign-flip trials to find the threshold for different ARL values. The results of WL-Sum, WL-Max, ST-Sum, and ST-Max procedures are presented in Table 3, indicating the approximation is reasonably accurate as the relative error is smaller than 1% for WL-Sum, around 7% to 11% for ST-Sum, around 1% for WL-Max and smaller than 8% for ST-Max. We also note that the magnitude of thresholds varies for different test statistics. For example, the WL-Sum statistic involves the summation of $p(p-1)/2$ terms and results in consistently higher thresholds than max-type statistics.

When a change point $\nu$ exists, the performance of the detection procedure is measured by the detection delay as defined in (2), which can be interpreted as the expected number of post-change samples needed to detect the change. We note that the supremum over all possible change-points $\nu$ in the definition (2) is typically intractable for a general detection procedure. Therefore, we use a simplified definition of the expected detection delay (EDD) as ${\mathbb{E}}_{1}[T]$, i.e., the expected stopping time when the change point equals $1$, this is also what is typically simulated in practice. We present the following proposition for the detection delay of the proposed stopping times.