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.
61,236 characters · 14 sections · 33 citation commands
Detection and Estimation of Structural Breaks in High-Dimensional Functional Time Series
\centerline{\bf Abstract}
In this paper, we consider detecting and estimating breaks in heterogeneous mean functions of high-dimensional functional time series which are allowed to be cross-sectionally correlated and temporally dependent. A new test statistic combining the functional CUSUM statistic and power enhancement component is proposed with asymptotic null distribution theory comparable to the conventional CUSUM theory derived for a single functional time series. In particular, the extra power enhancement component enlarges the region where the proposed test has power, and results in stable power performance when breaks are sparse in the alternative hypothesis. Furthermore, we impose a latent group structure on the subjects with heterogeneous break points and introduce an easy-to-implement clustering algorithm with an information criterion to consistently estimate the unknown group number and membership. The estimated group structure can subsequently improve the convergence property of the post-clustering break point estimate. Monte-Carlo simulation studies and empirical applications show that the proposed estimation and testing techniques have satisfactory performance in finite samples.
{\em Keywords}: clustering, CUSUM, functional time series, power enhancement, structural breaks.
\setcounter{equation}{0}
Modelling functional time series, time series of random functions defined within a finite interval, has became one of the main frontiers of developments in time series models. Various functional linear and nonlinear time series models have been proposed and extensively studied in the past two decades B00, HK10, HK12, HHR13, LRS20. These models together with relevant methodologies have been applied to various fields such as biology, demography, economics, environmental science and finance. However, the model frameworks and methodologies developed in the aforementioned literature heavily rely on the stationarity assumption, which is often rejected when testing the functional time series data in practice. For example, HKR14 find evidence of nonstationarity for intraday price curves of some stocks collected in the US market; ARS18 reject the null hypothesis of stationarity for the temperature curves collected in Australia; and LRS23 reveal evidence of nonstationary feature for the functional time series constructed from the age- and sex-specific life-table death counts. It thus becomes imperative to test whether the collected functional time series are stationary.
The primary interest of this paper is to test whether there exist structural breaks in the mean function over time and subsequently estimate locations of breaks if they do exist. There have been increasing interests on detecting and estimating structural breaks in functional time series. Broadly speaking, there are two types of detection techniques. One is to first reduce the infinite dimension of functional data to a finite dimension via the classic functional principal component analysis and then use the detection method developed for multivariate data to identify breaks AGHP09,BGHP09, ZSHW11, AK12, and the other is a fully functional detection method without preliminary dimension reduction HKR14, STW16, ARS18. The latter avoids possible information loss caused by the dimension reduction, thus may have more reliable numerical performance. The existing research works often focus on break detection in a single functional time series sequence. An extension to multiple functional time series is straightforward, but a further extension to high-dimensional functional time series is challenging. In fact, so far as we know, there is virtually no work on break detection and estimation in high-dimensional functional time series with dimension comparable to the sample size.
In this paper, we consider structural break detection within a flexible high-dimensional functional time series model framework, allowing the functional time series to be cross-sectionally correlated over subjects and temporally dependent over time. A novel test statistic combining the CUSUM statistic with the power enhancement (PE) component is proposed to detect breaks in heterogeneous mean functions. The asymptotic null distribution of the proposed test statistic is comparable to that of the conventional CUSUM statistic developed for a single functional time series HKR14, ARS18. In practice, it is often the case that breaks are sparse in the alternative hypothesis, i.e., breaks only occur in a small number of functional time series, or that the break locations may vary over subjects. Consequently, the traditional CUSUM test statistic would have low powers. To address this problem, an extra PE component is added to the test statistic. The PE technique is introduced by FLY15 to improve the power performance of high-dimensional tests for sparse alternatives and is further studied by KP19. We will show that incorporating the PE component in the break test enlarges the region where the proposed test has power, while avoids size distortion.
When a large number of functional time series are available, under the alternative hypothesis (with functional structural breaks), it is natural to split all the subjects into two groups: subjects of functional time series without breaks, and those with breaks. For the latter, we further estimate the unknown break points and derive a uniform approximation rate which is nearly optimal up to a logarithmic factor. In addition, we assume that there exists a latent group structure on the heterogeneous break points, i.e., subjects in each group have a common break point, whereas the break points are distinct between groups, and develop an easy-to-implement clustering method with an information criterion to consistently estimate the latent structure including the group membership and number of distinct break points. This complements the recent literature on cluster analysis of functional data which are often assumed to be independent and identically distributed (i.i.d.) or stationary and weakly dependent TK03, FV06, CL07, PM08, DH12, DHP19. Furthermore, we estimate the homogenous break point within each group and show that the consistency property can be improved for the estimated break point by making use of the estimated group structure.
The Monte-Carlo simulation studies show that the proposed test with the extra PE component achieves the power enhancement when structural breaks are relatively sparse while avoids size distortion in finite samples. We also compare the numerical performance between different choices of high-criticism thresholding parameters. With the latent group structure on change points, the proposed clustering algorithm and information criterion can accurately identify the cluster number and membership. Furthermore, the pooled CUSUM method using the estimated cluster structure can substantially improve the estimation accuracy of the common break points. The developed break detection and estimation methods are applied to Dow Jones Industrial Average constituent stocks and the age-specific mortality rates collected in $32$ countries to test the existence of structural breaks and estimate the break locations. The empirical analysis reveals that there exists one common break for either of the two real data sets.
The rest of the paper is organised as follows. Section (ref) presents the model framework, hypothesis testing problem and some fundamental assumptions. Section (ref) introduces the test statistic together with its asymptotic properties under both the null and alternative hypotheses. Section (ref) considers the break point estimation and identifies the latent structure on the heterogeneous break points. Sections (ref) and (ref) report the simulation and empirical studies, respectively. Section (ref) concludes the paper. Proofs of the main asymptotic theorems are given in Appendix A whereas proofs of some technical lemmas and propositions are available in a supplement. Throughout the paper, we let ${\mathscr H}$ be the Hilbert space defined as a set of measurable functions $f(\cdot)$ on a bounded set ${\mathbb C}$ such that $\int_{\mathbb C} f^2(u)du<\infty$. The inner product in ${\mathscr H}$ is defined as $\langle f_1,f_2\rangle=\int_{\mathbb C}f_1(u)f_2(u)du$. For $f\in{\mathscr H}$, we define the $L_2$-norm $\Vert f\Vert:=\Vert f\Vert_2=\langle f,f\rangle^{1/2}=\left[\int_{\mathbb C}|f(u)|^2du\right]^{1/2}$ and more generally the $L_p$-norm $\Vert f\Vert_p=\left[\int_{\mathbb C}|f(u)|^pdu\right]^{1/p}$, $p\geq1$. Let ${\mathscr L}({\mathscr H})$ be the space of continuous linear operators from ${\mathscr H}$ to ${\mathscr H}$ equipped with the operation norm defined by $\Vert{\mathcal L}\Vert_O=\sup_{f\in{\mathscr H}}\left\{\Vert{\mathcal L}(f)\Vert: \Vert f\Vert\leq1\right\}$. For each ${\mathcal L}\in{\mathscr L}({\mathscr H})$, its adjoint ${\mathcal L}^\prime$ is defined via $\langle {\mathcal L}f_1,f_2\rangle=\langle f_1,{\mathcal L}^\prime f_2\rangle$ for any $f_1,f_2\in{\mathscr H}$. Let $\stackrel{d}\longrightarrow$, $\stackrel{P}\longrightarrow$ and $\Rightarrow$ denote convergence in distribution, convergence in probability and weak convergence, respectively.
\setcounter{equation}{0}
Suppose that we collect a sequence of functional observations ${\mathbf X}_t=(X_{1t},\cdots,X_{Nt})^{^\intercal}$, $t=1,\cdots,T$, where $X_{it}=\left(X_{it}(u): u\in{\mathbb C}\right)$. For the $i$-th subject, $X_{it}$, $t=1,\cdots,T$, are generated from the following model with a possible break in the mean function:
where $\mu_i=(\mu_i(u): u\in{\mathbb C})$ is the pre-break mean function, $\delta_i=(\delta_i(u): u\in{\mathbb C})$ is the jump function, $\tau_i$ is the break point, $I(\cdot)$ is the indicator function and $\epsilon_{it}=(\epsilon_{it}(u): u\in{\mathbb C})$ is stationary over the temporal dimension. The unknown functional components $\mu_i$ and $\delta_i$ as well as the break points $\tau_i$ are allowed to vary over $i$, reflecting heterogeneity of functional time series over subjects. Our primary interest lies in the following hypothesis testing problem:
Under the null hypothesis $H_0$, model (ref) reduces to $X_{it}=\mu_i+\epsilon_{it}$, which is stationary over $t$. We aim to test whether to reject $H_0$ and estimate the break points $\tau_i$ if $H_0$ is rejected. The functional errors $\epsilon_{it}$ satisfy the following regularity conditions.
\setcounter{assumption}{0}
From Assumption (ref)(i), the functional linear process $\epsilon_{it}$ defined in (ref) with coefficient operators satisfying (ref) is stationary and short-range dependent over time $t$, but its distribution is allowed to vary over $i$. This assumption is a natural extension of the assumptions in Ba10 and HH12 for panel time series setting. Assumption (ref)(ii) shows that the functional time series may be weakly cross-sectionally correlated over $i$. The moment restriction in (ref) is a high-level condition which is easy to verify when $\eta_{it}$ are independent over $i$. In the latter case, we may further show that \[\widetilde{\boldsymbol\Omega}=\lim_{N\rightarrow\infty}\frac{1}{N}\sum_{i=1}^N{\mathbf A}_i{\boldsymbol\Omega}_i{\mathbf A}_i^\prime,\] where ${\boldsymbol\Omega}_i$ is the integral operator with $\Omega_i(u,v)$ as the kernel. As $\widetilde{\boldsymbol\Omega}$ is positive definite, we may conduct an eigenanalysis and find pairs of non-negative eigenvalues and eigenvectors $(\lambda_{k},\psi_{k})$ (with eigenvalues arranged in an non-increasing order), $k=1,2,\cdots$, such that
The exponential moment condition ${\sf E}\left[\exp\left\{c_\eta\Vert\eta_{it}\Vert^2\right\}\right]<\infty$ in Assumption (ref)(ii) facilitates the truncation technique and application of the concentration inequality for random elements in ${\mathscr H}$ B00, BLM16 in proofs of the main asymptotic results.
\setcounter{equation}{0}
Letting $\nu_{it}=\mu_i+\delta_i I\left(t>\tau_i\right)$, we may re-write model (ref) as
Define $$\widetilde{X}_t=\frac{1}{N}\sum_{i=1}^N X_{it},\ \ \widetilde{\nu}_t=\frac{1}{N}\sum_{i=1}^N\nu_i\ \ {\rm and}\ \ \widetilde{\epsilon}_t=\frac{1}{N}\sum_{i=1}^N \epsilon_{it},$$ where we suppress their dependence on $N$. From the model formulation (ref), we have
Note that, under $H_0$, the time-varying mean function $\widetilde{\nu}_t$ becomes $\widetilde{\mu}=\frac{1}{N}\sum_{i=1}^N\mu_{i}$ which is time-invariant and model (ref) thus reduces to \[\widetilde{X}_t=\widetilde{\mu}+\widetilde{\epsilon}_t,\ \ t=1,\cdots,T.\] In order to test $H_0$ in (ref), a naive idea is to construct the CUSUM test statistic using $\widetilde{X}_{t}$. Define the functional CUSUM statistic as
where $0\leq x\leq 1$, $u\in{\mathbb C}$ and $\lfloor \cdot\rfloor$ denotes the floor function, and subsequently construct the test statistic via
The $\sqrt{N}$-rate in the definition of $\widetilde{Z}_{NT}(x;u)$ is appropriate when $X_{it}$ and $\epsilon_{it}$ are weakly dependent (or independent as a special case) over $i$. If the functional time series are strongly dependent over $i$, a different normalisation rate may be required. Under some regularity conditions, we may show that $Z_{NT}$ has the same asymptotic null distribution as the classic CUSUM test statistic for a single functional time series, see Theorem (ref) below. However, the test statistic based on the model formulation (ref) suffers from the low-power issue as $\widetilde{\nu}_t$ may be time-invariant even when there are significant breaks in the subject-specific mean functions $\nu_{it}$ for some indices $i$. For example, in the case of common break ($\tau_i\equiv \tau$), when $N$ is even, $\delta_i=\delta$ for $i=1,\cdots,N/2$ and $\delta_i=-\delta$ for $i=N/2+1,\cdots,N$, we have $\widetilde\nu_t=\widetilde\mu$. Consequently, $\widetilde{X}_{t}=\tilde{\mu}+\widetilde{\epsilon}_{t}$ under both $H_0$ and $H_A$ and the CUSUM test statistic in (ref) would have no power. If $\delta_i\neq 0$ only holds for a fixed number of indices $i$ (i.e., breaks are sparse), we may show that $\widetilde\nu_t\approx\widetilde\mu$, which would also result in low power of the CUSUM test.
To fix the low-power issue of the functional CUSUM test statistic defined in (ref) and (ref), we next consider incorporating a PE component in construction of the break test statistic. FLY15 introduce the PE technique in high-dimensional cross-sectional tests, improving the power performance in testing sparse-type alternatives. The PE component is expected to satisfy the following three properties: (i) non-negativity, (ii) no size-distortion, and (iii) power enhancement. We aim to augment the CUSUM test statistic by adding a PE component which satisfies these three properties. Define the subject-specific functional CUSUM statistic:
and the PE component:
where $\xi_{NT}$ denotes a high-criticism threshold and $a\vee b=\max\{a,b\}$. The PE component $Z_{NT}^\diamond$ is clearly non-negative. We need to show that ${\sf P}\left(Z_{NT}^\diamond=0 | H_0\right)\rightarrow1$ and $Z_{NT}^\diamond$ diverges to infinity with probability approaching one ({\em w.p.a.1}) under some specific regions of the alternative hypothesis on which the CUSUM test defined in (ref) and (ref) has low power, verifying properties (ii) and (iii) of the PE component (see the proof of Theorem (ref)). In particular, we adopt the diverging factor $\sqrt{N\vee T}$ in the PE component rather than $\sqrt{N}$ used by FLY15, to ensure that property (iii) can be achieved no matter $N$ is fixed or diverging. The formulation of $Z_{NT}^\diamond$ in (ref) is analogous to the sparsified CUSUM statistic proposed by CF15 for high-dimensional real-valued time series. Throughout this paper we set the threshold as $\xi_{NT}=c_\xi\ln(N\vee T)\ln\ln (N\vee T)$ with $c_\xi$ being a user-specified positive constant. In Section (ref) below, we also consider $\xi_{NT}=c_\xi\ln(NT)\ln\ln (NT)$ in the simulation study and compare the finite-sample numerical performance between these two choices of high-criticism thresholding parameters.
Combining the CUSUM test statistic $Z_{NT}$ defined in (ref) and the PE component $Z_{NT}^\diamond$ defined in (ref), we propose the following power enhanced CUSUM (PE-CUSUM) test statistic:
Let $H_A^\diamond$ be the alternative $H_A$ defined in (ref) with $\delta_i$ satisfying
and $\widetilde{H}_A^\diamond$ the alternative $H_A$ with $\delta_i$ satisfying (ref) or $\widetilde{\nu}_t$ defined in (ref) satisfying
The following theorem gives the asymptotic properties of the PE-CUSUM test statistic $\widehat{Z}_{NT}$.
\setcounter{theorem}{0}
The condition $N=O(T^{\kappa})$ indicates that the number of subjects $N$ can be much larger than the time series length $T$ ($\kappa>1$). The asymptotic null distribution of the PE-CUSUM test statistic is similar to that in Theorem 2.1 of HKR14, Corollary 1 of STW16 and Theorem 1 of ARS18. The involvement of the PE component does not lead to size distortion, i.e., ${\sf P}\left(Z_{NT}^\diamond=0 | H_0\right)\rightarrow1$ and property (ii) of the PE component is thus satisfied. The limit distribution in (ref) relies on the unknown eigenvalues $\lambda_i$ and the standard Brownian bridges $B_i(x)$. STW16 introduce a block bootstrap method to choose the critical value of the test, whereas ARS18 suggest a Monte-Carlo simulation method.
Theorem (ref)(ii) shows that the developed test is consistent under the alternative hypothesis when either (ref) or (ref) is satisfied. Without the PE component $Z_{NT}^\diamond$, the CUSUM test statistic is only consistent under the restrictive high-level condition (ref), which can be seen as a diverging time-varying measure of the mean functions. Note that the latter condition is often violated when breaks from different subjects are cancelled out as discussed earlier in the section. By incorporating $Z_{NT}^\diamond$, the PE property is achieved, i.e., the region where the proposed test has power is enlarged from $H_A$ with (ref) to $H_A$ with either (ref) or (ref). In particular, the condition (ref) covers the scenario of sparse breaks. It follows from (ref) that the PE property is achieved when $N$ is either fixed or diverging.
\setcounter{equation}{0}
When $H_0$ is rejected, it is often the case that structural breaks only occur in some of the functional time series processes. Hence, we may split the index set $\{1,2,\cdots,N\}$ into \[{\mathcal C}_\bullet=\{ 1\leq i\leq N: \delta_i\neq 0\}\ \ {\rm and}\ \ {\mathcal C}_\circ=\{ 1\leq i\leq N: \delta_i=0\},\] the index set with structural breaks in mean functions and that without breaks. It is natural to estimate ${\mathcal C}_\bullet$ and ${\mathcal C}_\circ$ by \[ \widehat{\mathcal C}_\bullet=\left\{ 1\leq i\leq N: \sup_{0\leq x\leq 1}\int_{\mathbb C} Z_{iT}^2(x;u)du\geq\xi_{NT}\right\}\] and \[\widehat{\mathcal C}_\circ=\left\{ 1\leq i\leq N: \sup_{0\leq x\leq 1}\int_{\mathbb C} Z_{iT}^2(x;u)du<\xi_{NT}\right\}, \] respectively, where $Z_{iT}(x;u)$ and $\xi_{NT}$ are defined in Section (ref). Subsequently, we estimate the (heterogeneous) break points $\tau_i$ by
The following theorem shows that $\widehat{\mathcal C}_\bullet$ and $\widehat{\mathcal C}_\circ$ are consistent estimates and provides a uniform approximation rate for $\widehat{\tau}_i$ over $i\in{\mathcal C}_\bullet$.
As $N=O(T^\kappa)$, the uniform approximation rate in (ref) can be simplified to $(\ln T)^{1+\zeta}$. The uniform approximation rate for the estimated break point is nearly optimal up to a logarithmic factor, and is slightly slower than the $O_P(1)$ rate derived for the single functional time series setting. We next show that the approximation rate in (ref) can be further improved by imposing a latent group structure on the break points $\tau_i$, $i\in{\mathcal C}_\bullet$.
Assume that there exists a partition of ${\mathcal C}_\bullet$, denoted by ${\mathcal C}(b_1),\cdots,{\mathcal C}(b_{K_0})$, such that
where $b_1<b_2<\cdots<b_{K_0}$ are $K_0$ distinct break points. Neither the group membership nor the number $K_0$ is known. Note that the common break setting is nested in the latent structure (ref) with $b_1=\cdots=b_{K_0}=\tau_0$ and ${\mathcal C}(\tau_0)={\mathcal C}_\bullet$. The latent structure (ref) shows that $X_{it}$ have the common break point $b_k$ when $i\in{\mathcal C}(b_k)$.
With the estimated break points $\widehat{\tau}_i$ defined in (ref), $i\in\widehat{\mathcal C}_\bullet$, we first sort them from minimum to maximum and denote the ordered points as $\widehat{\tau}^{(1)},\cdots,\widehat\tau^{(n)}$, $n=\big|\widehat{\mathcal C}_\bullet \big|$, where $|{\cal A}|$ denotes the cardinality of set ${\cal A}$. Then calculate the jumps: \[ \Delta_i(\widehat{\tau})=\widehat{\tau}^{(i+1)}-\widehat{\tau}^{(i)},\ \ i=1,\cdots,n-1. \] If the number of distinct break points is assumed to be $K$, define $\overline{\tau}_j=\widehat{\tau}^{(i+1)}$ with $\Delta_i(\widehat{\tau})$ being the $j$-th largest jump, $1\leq j\leq K-1$. Then we sort $\overline{\tau}_1,\cdots,\overline{\tau}_{K-1}$ from minimum to maximum, denote them as $\overline{\tau}_K^{(1)},\cdots,\overline{\tau}_K^{(K-1)}$, and obtain the estimated clusters as
where, without loss of generality, $\overline{\tau}_K^{(0)}=1$ and $\overline{\tau}_K^{(K)}=T$. Given the cluster number $K$, we compute
with \[ \widehat{\mu}_{i,k|K}=\frac{1}{\widehat{\tau}_{k|K}}\sum_{t=1}^{\widehat{\tau}_{k|K}}X_{it}\ \ \ {\rm and}\ \ \ \widehat{\delta}_{i,k|K}=\frac{1}{T-\widehat{\tau}_{k|K}}\sum_{t=\widehat{\tau}_{k|K}+1}^TX_{it}-\frac{1}{\widehat{\tau}_{k|K}}\sum_{t=1}^{\widehat{\tau}_{k|K}}X_{it}. \] Construct the penalised objective function:
where $\rho_{NT}$ is a user-specified tuning parameter satisfying some restrictions, see Assumption (ref)(iii) below, and \[ V(K)= \frac{1}{\big|\widehat{\mathcal C}_\bullet\big|}\sum_{k=1}^K\sum_{i\in\widehat{\mathcal C}(k|K)}\frac{1}{T}\sum_{t=1}^T\left\Vert X_{it}-\widehat{\nu}_{it,k|K}\right\Vert^2. \] The true cluster number $K_0$ is determined by
where $\overline{K}$ is a pre-specified upper bound of the cluster number. Replacing $K$ by $\widehat{K}$ in (ref), we obtain \[ \widehat{\mathcal C}(b_k)=\widehat{\mathcal C}\big(k|\widehat{K}\big),\quad k=1,\cdots,\widehat{K}. \] The following assumption is required to derive the consistency property for $\widehat{K}$ and $\widehat{\mathcal C}(b_k)$.
Assumption (ref)(i) indicates that the distance between distinct break points is of the same order as $T$ and the cardinality of ${\mathcal C}(b_k)$ is similar over $k=1,\cdots,K_0$. Assumption (ref)(ii) imposes some mild restrictions on sizes of $\Vert\mu_i\Vert$ and $\Vert\delta_i\Vert$. Finally, Assumption (ref)(iii) is a crucial condition on the tuning parameter in the penalty term, ensuring that the information criterion can consistently select $K_0$.
The above theorem establishes the consistency property for the proposed information criterion and cluster analysis of the heterogeneous break points. As $X_{it}$, $i\in{\mathcal C}(b_k)$, $t=1,\cdots,T$, have a common break point, it is sensible to estimate $b_k$ more efficiently by pooling the CUSUM quantities over the subjects in $\widehat{\mathcal C}(b_k)$, i.e.,
whose consistency property is given in the following theorem.
Theorem (ref) extends Theorem 3.1 in Ba10 to the high-dimensional functional data setting with a latent structure on the heterogeneous break points. The condition (ref) plays a key role in the theoretical derivation and is a natural extension of Assumption 2 in Ba10. It indicates that break sizes cannot be too small over $i\in{\mathcal C}_\bullet$ so that the true break point can be consistently estimated. In particular, when $\min_{i\in{\mathcal C}_\bullet}\Vert\delta_i\Vert$ is bounded away from zero as in Assumption (ref)(ii), (ref) would be automatically satisfied if $|{\mathcal C}_\bullet|\rightarrow\infty$. Note that in the case of break estimation for a single time series process, the optimal asymptotic order for the break point estimation is $O_P(1)$ Ba97, ARS18. With more sample information from pooling a panel of functional time series, it is unsurprising to improve the approximation order to $o_P(1)$. The restrictions $|{\mathcal C}_\bullet |=O(T^2)$ and $T=O(|{\mathcal C}_\bullet|^{3/2})$ indicate that $|{\mathcal C}_\bullet|$ may diverge to infinity at a faster rate than $T$. In addition, the cross-sectional independence condition on $\eta_{it}$ facilitates the technical proofs but may be replaced by some high-level conditions such as those in Assumption (ref)(ii).
\setcounter{equation}{0}
Generate time series of random functions $[\epsilon_{i1}(u),\epsilon_{i2}(u),\cdots,\epsilon_{iT}(u)]$, $i=1,\cdots,N$, as follows
where $f_1(u)$, $f_2(u), \cdots, f_{J}(u)$ are randomly sampled without replacement from $J=21$ Fourier basis functions, and $\eta_{it,j}$ are innovations following ${\sf N}(0,1/j)$ independently over $i$, $t$ and $j$. To generate functional values, we consider $101$ equally-spaced grid points between $0$ and $1$. Writing $\boldsymbol{\beta}_{t,j} = \left(\beta_{1t,j},\beta_{2t,j},\cdots,\beta_{Nt,j}\right)^{^\intercal}$, $j=1,\cdots,21$, we generate $\boldsymbol{\beta}_{t,j}$ (independently over $j$) from a vector autoregression of order 1:
where $\boldsymbol{A}=(a_{ij})_{N\times N}$ is the transition matrix, and $\boldsymbol{z}_{t,j}$ is independently generated by an $N$-dimensional normal distribution with mean zero and identity covariance matrix. Following LRS23, $\boldsymbol{A}$ is a banded matrix with $a_{ij}$ independently generated from a ${\sf U}(-0.3,0.3)$ when $|i-j|\leq 3$ and $a_{ij}=0$ when $|i-j|>3$. In the simulation, we set $T=200$, $N=200$ or $400$, and the replication number $R=1000$.
We randomly select ${\sf SDR} \times N$ subjects with a change point, where {\sf SDR} denotes a sparse-to-dense ratio. The remaining subjects have no change point. To specify a change-point location $\tau_i$ for the selected subject, we draw ${\sf SDR}\times N$ values from ${\sf U}(0.25 \times T, 0.75\times T)$. The lower and upper bounds of the uniform distribution are purposely chosen so that the change point location is not close to the boundary of a sample. As in ARS18, for $i\in{\mathcal C}_\bullet$ (the index set for subjects with breaks), we define a class of break functions: \[ \delta_{i,m}(u) = \delta_m^*(u)\times \sqrt{c_{i}^\ast},\ \ \delta_m^*(u)= \frac{1}{\sqrt{m}}\sum^m_{j=1}f_{j}(u), \] where $m=1,\cdots,J$ and $c_i^\ast$ is a positive constant to be specified later. For each $i\in{\mathcal C}_\bullet$, $\delta_{i,1}(u)$ is the case of a break only in the leading eigendirection (determined by the basis function $f_1$), while $\delta_{i,J}(u)$ is the case of a break that affects all the eigendirections, see the discussion in ARS18. The value $c_i^\ast$ controls the break magnitude, linking to the signal-to-noise ratio: \[ \text{SNR}_{i} = \frac{(\tau_i/T)(1-\tau_i/T)\|\delta_{i,m}\|^2}{\text{tr}({\boldsymbol\Omega}_{\epsilon_{i}})}=c_{i}^\ast\times\frac{(\tau_i/T)(1-\tau_i/T)}{\text{tr}({\boldsymbol\Omega}_{\epsilon_{i}})}, \] where $\text{tr}(\cdot)$ denotes the trace of a square matrix and ${\boldsymbol\Omega}_{\epsilon_{i}}$ denotes the long-run covariance matrix of ${\boldsymbol\epsilon}_{it}=\left(\epsilon_{it,1},\epsilon_{it,2},\cdots, \epsilon_{it,21}\right)^{^\intercal}$ over the time span. The value of $c_{i}^\ast$ can be easily computed with a given SNR$_{i}$ level. For those subjects without breaks, we set $c_{i}^\ast=0$. With a given integer value of $m$, we finally simulate a panel of functional time series as follows,
In the simulation study, we consider the following two choices of high-criticism thresholds in the PE-CUSUM test statistic:
where $c_\xi=\lambda_1^{1/2}$ with $\lambda_1$ being the leading eigenvalue of $\widetilde{\boldsymbol\Omega}$ in (ref). The asymptotic theorems in Sections (ref) and (ref) are derived by setting $\xi_{NT}=\xi_{NT,2}$. In fact, they continue to hold with minor modifications when $\xi_{NT}=\xi_{NT,1}$. In the following simulation results, we denote these two PE-CUSUM test statistics as PE-CUSUM$_1$ and PE-CUSUM$_2$. We also consider the conventional CUSUM (without PE) defined in (ref) as a benchmark. We only report the simulation results when $m=1$ in the break function definition, i.e., breaks occur in the leading eigendirection, since the results are similar when $m$ is set as other positive integer values.
In Table (ref) below, we report the size performance of CUSUM and PE-CUSUM in finite samples, where three levels of significance $\alpha=0.01, 0.05$, and 0.10 are considered. The test critical values are determined by ARS18's simulation-based method. The sizes of all the three tests are generally close to the nominal ones. In particular, the results confirm the validity of Theorem (ref)(i) and show that incorporating the PE component in the CUSUM test statistic does not lead to severe size distortion.
In Table (ref), we report the power performance of the CUSUM and PE-CUSUM test statistics when SNRs are $10^{-1}, 10^{-2}$ and $5\times 10^{-3}$ and SDRs are $0.1$ and $0.5$. We note that percentages of rejecting the null hypothesis via the three test statistics are the same and close to $1$ when SDR is $0.5$, whereas both the PE-CUSUM$_1$ and PE-CUSUM$_2$ outperform the CUSUM statistic (without PE) when breaks are sparse (i.e., SDR is $0.1$) and the signal-to-noise ratios are low (i.e., SNRs are $10^{-2}$ and $5\times 10^{-3}$). The latter confirms the improvement of power performance with the extra PE component in the test statistic (ref). Meanwhile, the performance of PE-CUSUM$_2$ is more stable than PE-CUSUM$_1$ in particular when SNR is $5\times 10^{-3}$, indicating that $\xi_{NT,2}= c_\xi\ln(N\vee T) \ln\ln(N\vee T)$ may be a more appropriate choice for the high-criticism threshold in PE-CUSUM.
{
}
We next compute the percentage of correctly identifying subjects with a change point, i.e, ${\mathcal C}_\bullet$ defined in Section (ref). It is defined by $\text{TP}/N$, where TP denotes the number of true positive hits (the subject indices match with the estimated indices), and $N$ denotes the number of subjects. Let the subjects with breaks be categorised as one, and those without breaks be categorised as zero. We also compute the $F_1$ score, which is the harmonic mean of precision defined by \[ \text{F}_1 = \frac{\text{TP}}{\text{TP}+(\text{FP}+\text{FN})/2}, \] where FP denotes false alarm or overestimation (type I error, where a change point is detected for subjects without breaks), and FN denotes miss or underestimation (type II error, where a change point is not detected for subjects with breaks).
The TP and $\text{F}_1$ score results are reported in Table (ref), when SNRs are $10^{-1}, 10^{-2}$ and $5\times 10^{-3}$ and SDRs are $0.1$ and $0.5$. As SNR decreases from $10^{-1}$ to $5\times 10^{-3}$, there is a decrease of the TP and F$_1$ values. This is unsurprising as it becomes more difficult to detect breaks when the signal-to-noise ratio is lower. The estimation $\widehat{\mathcal C}_\bullet$ using $\xi_{NT}=\xi_{NT,2}$ performs better than that with $\xi_{NT}=\xi_{NT,1}$, which is consistent with the discussion in Section (ref), and again justifies the use of $\xi_{NT}= c_\xi\ln(N\vee T) \ln\ln(N\vee T)$ as the high-criticism threshold.
We next assess the clustering methodology proposed in Section (ref) by imposing a latent group structure on change points. Specifically, we set the number of distinct change points as $K_0=3$ and split the $N$ subjects into the following four groups: the first half of $N$ subjects have no change point; and the remaining subjects are equally split into three groups with the change point at $0.25\times T$, $0.5\times T$ and $0.75\times T$, respectively. The tuning parameter used in the information criterion is set as $\rho_{NT}=(N\vee T)^{-1/2} \ln(N \vee T)$, so that Assumption (ref)(iii) is satisfied (if $N=T^\phi$ with $\phi<2$). To assess the estimation accuracy of group membership, we compute the Purity and normalised mutual information (NMI) measurements which are respectively defined as \[ \text{Purity}\left(\widehat{\mathcal C}_\bullet, {\mathcal C}_\bullet\right) = \frac{1}{N}\sum^{\widehat{K}}_{k=1}\max_{1\leq j\leq K_0}\left|\widehat{\mathcal{\mathcal C}}_k\cap {\mathcal C}_j\right|, \] and \[ \text{NMI}\left(\widehat{\mathcal C}_\bullet, {\mathcal C}_\bullet\right) = 2\frac{ I\left(\widehat{C}_\bullet, {\mathcal C}_\bullet\right)}{H(\widehat{\mathcal C}_\bullet)+H({\mathcal C}_\bullet)}, \] where $\widehat{\mathcal C}_\bullet = \left\{\mathcal{\mathcal C}_1,\cdots,\mathcal{\mathcal C}_{\widehat{K}}\right\}$ is the estimate of ${\mathcal C}_\bullet=\left\{\mathcal{C}_1,\cdots, \mathcal{C}_{K_0}\right\}$, $H({\mathcal C}_\bullet)$ denotes the entropy of ${\mathcal C}_\bullet$, $I\left(\widehat{\mathcal C}_\bullet, {\mathcal C}_\bullet\right)$ is the mutual information between $\widehat{\mathcal C}_\bullet$ and ${\mathcal C}_\bullet$ defined by \[ I\left(\widehat{\mathcal C}_\bullet, {\mathcal C}_\bullet\right) = \sum^{\widehat{K}}_{k=1}\sum^{K_0}_{j=1}\left(\frac{|\widehat{\mathcal{C}}_k \cap \mathcal{C}_j|}{N}\right)\log_2\left(\frac{N|\widehat{\mathcal{C}}_k\cap {\mathbb C}_j|}{|\widehat{\mathcal{C}}_k||\mathcal{C}_j|}\right). \]
The relevant results are summarised in Tables (ref) and (ref). It follows from Table (ref) that the frequency of correctly estimating the cluster number generally decreases as the SNR decreases from $10^{-1}$ to $5\times 10^{-3}$ (because it becomes more difficult to detect breaks when SNR is lower). The information criterion with $\xi_{NT}=\xi_{NT,2}$ can more accurately estimate the cluster number than that with $\xi_{NT}=\xi_{NT,1}$. In particular, when SNRs are $10^{-2}$ and $5\times 10^{-3}$, the information criterion with $\xi_{NT}=\xi_{NT,1}$ tends to under-estimate the cluster number. The Purity and NMI values in Table (ref) are close to the perfect value of one when SNR is $10^{-1}$, and gradually decrease when SNR becomes smaller.
We finally assess the performance of the estimated break locations, and compare the post-clustering estimation (ref) with the pre-clustering estimation (ref), from which we may demonstrate the usefulness of clustering heterogeneous change points over subjects. We compute the mean squared distances between the estimated and true change points with the results reported in Table (ref). According to the previous simulation results, it may be more appropriate to use $\xi_{NT}=\xi_{NT,2}$ in the break location estimation. Table (ref) shows that, conditional on the accurate estimation of the cluster number, the post-clustering estimation performs significantly better than the pre-clustering one which ignores the latent group structure on the heterogeneous change points. This is consistent with the convergence results in Theorems (ref) and (ref), which state that the post-clustering estimation is consistent, whereas the pre-clustering estimation has the logarithmic-($N\vee T$) approximation order.
\setcounter{equation}{0}
In this section, we apply the developed break detection and estimation methods to two empirical data sets: $28$ Dow Jones Industrial Average (DJIA) constituent stocks from January 2, 2018 to December 31, 2021, and the age-specific mortality rates collected in $32$ countries from 1960 to 2013.
The DJIA index shows how $30$ publicly owned large companies based in the United States have traded during a standard New York Stock Exchange trading session. Table (ref) lists the stock names and tick symbols of 30 constituents for the DJIA index. We consider their daily cross-sectional returns from January 2, 2018 to December 31, 2021, with the data obtained from the Refinitiv Datascope (\url{https://select.datascope.refinitiv.com/DataScope/}). There are $T=1,008$ trading days. Among the 30 constituent stocks, DOW and HON.O began trading on April 2, 2019 and May 11, 2021, respectively. Thus, we remove these two stocks in our empirical analysis, resulting in $N=28$. For each trading day, we consider $5$-minute resolution data covering the period between 9:30 and 15:55 Eastern standard time, and obtain $78$ data points. For asset $i$, let $P_{it}(u_j)$ be the intraday 5-minute close price at time $u_j$ on trading day $t$, and construct a sequence of CIDRs RWZ20: \[ X_{it}(u_j) = 100\times [\ln P_{it}(u_j) - \ln P_{it}(u_1)], \qquad j=2,3,\cdots,78, \] where $i=1,2,\cdots,28$ and $t=1,\cdots,1,008$. We use the linear interpolation algorithm Hyndman19 to convert discrete data points into a continuous function.
We implement the proposed PE-CUSUM test to detect if there exists a structural break in at least one of the $28$ constituents. The test $p$-value is 0, indicating the existence of breaks. We then estimate the heterogeneous change points using (ref) and identify the constituents with a structural break. The estimation results are reported in Table (ref). From the estimated break dates, it is sensible to expect a common change point (likely in March 2023) shared by some stocks. Hence, we further implement the clustering algorithm proposed in Section (ref) to estimate the number of clusters and common break date. The information criterion (ref) selects the number of clusters as one, and the post-clustering pooled CUSUM method in (ref) estimates the common break date on March 16, 2020. This is consistent with our observation of the heterogeneous break date estimates in Table (ref). The common break is related to the stock market crash in March 2020 when the DJIA index suffered severe losses on March 9 (-7.79%), March 12 (-9.99%), and March 16 (-12.93%).
The age-specific mortality rates are obtained from HMD23. Our dataset covers the period from 1960 to 2013 and $32$ countries with sufficient data to use. Table (ref) shows a list of these countries and the corresponding ISO Alpha-3 codes. For each of the 32 countries, we consider the ages from 0 to 99 in a calendar year and the last age group 100+. We smooth the age-specific mortality rates by a weighted penalised regression with monotonic constraint Wood94, HU07. The same data set is also considered by TSY22.
The developed PE-CUSUM test rejects the null hypothesis of no structural break. For the female population, Japan is the only country with a structural break and the estimated change point is 1983. For the male population, we use (ref) to estimate the heterogeneous change points and identify the countries with a structural break. The estimation results are reported in Table (ref), where we find that many developed countries experienced a structural break. It seems that a common break may occur in the late 1980s for those countries. This is confirmed by implementing the proposed clustering algorithm: the information criterion determines one cluster and the post-clustering CUSUM estimate of the common change point is 1988.
We propose a new fully functional test statistic combining the classic CUSUM and an extra PE component to detect structural breaks in the heterogeneous mean functions for large-scale functional time series, where the number of subjects may be larger than the time series length. The underlying functional time series are allowed to be weakly correlated over subjects. We derive the asymptotic property for the developed test under both the null and alternative hypotheses. In particular, the involvement of the PE component in the test statistic can enlarge the region where the test has power and detect sparse breaks in the alternative. We further impose a latent group structure on the heterogeneous break points and combine a simple clustering algorithm with an information criterion to accurately estimate the group membership and number. The post-clustering pooled CUSUM method using the estimated group structure is introduced to consistently locate the homogenous break point within each group. The Monte-Carlo simulation studies demonstrate the power enhancement property of the developed test in finite samples when breaks are sparse, the accuracy of the latent group structure estimation, and the convergence improvement of the post-clustering break point estimation (over the estimation neglecting the latent structure). The developed methodology is applied to detect and estimate structural breaks for DJIA constituent stocks and the age-specific mortality rates collected in $32$ countries.
The first author is partially supported by the Australian Research Council Discovery Project (DP230102250) and the National Natural Science Foundation of China (72033002). The second author is partially supported by the National Science Foundation grant (DMS 1820702) and National Institutes of Health grants (R01AI136664 and R01AI170249). The third author is partially supported by the Australian Research Council Discovery Project (DP230102250). The usual disclaimer applies.