{'fraction_non_alphanumeric': 0.05929183653815287, 'fraction_numerical': 0.021058714005970917, 'mean_word_length': 3.9380201315685186, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 37, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
High-dimensional time series segmentation via factor-adjusted vector autoregressive modelling

Vector autoregressive (VAR) models are popularly adopted for modelling high-dimensional time series, and their piecewise extensions allow for structural changes in the data. In VAR modelling, the number of parameters grow quadratically with the dimensionality which necessitates the assumption of sparsity in high dimensions. However, it is debatable whether sparse VAR models are adequate for handling datasets exhibiting strong serial and cross-sectional correlations. We propose a piecewise stationary time series model that simultaneously allows for strong correlations as well as structural changes, where pervasive serial and cross-sectional correlations are accounted for by a (possibly) time-varying factor structure, and any remaining idiosyncratic dependence between the variables is handled by a piecewise stationary, sparse VAR model. We propose an accompanying two-stage change point detection methodology which fully addresses the challenges arising from not observing either the factor-driven or the VAR processes directly. Its consistency in estimating both the total number and the locations of the change points in the latent components, is established under conditions considerably more general than those in the existing literature. We demonstrate the competitive performance of the proposed methodology on simulated datasets and an application to US blue chip stocks data.


Introduction
Vector autoregressive (VAR) models have emerged as a popular approach to model crosssectional and serial correlations in multivariate, possibly high-dimensional time series. See, for example, recent contributions in finance (Barigozzi and Hallin, 2017; Barigozzi and Brownlees, 2019), biology (Shojaie and Michailidis, 2010) and genomics (Michailidis and d'Alché Buc, 2013). Within such settings, the importance of data segmentation is wellrecognised, and several methods exist for detecting change points under piecewise stationary VAR models in both fixed (Kirch et al., 2015) and high dimensions (Safikhani and Shojaie, 2022;Wang et al., 2019;Bai et al., 2020;Maeng et al., 2022). VAR modelling quickly becomes a high-dimensional problem as the number of parameters grows quadratically with the dimensionality, and the assumption of sparsity is commonly made for estimation and inference (Basu and Michailidis, 2015;Davis et al., 2016;Liu and Zhang, 2021). Accordingly, most existing methods for detecting change points in high-dimensional, piecewise stationary VAR processes require the sparsity of VAR parameters. However, it is debatable whether highly sparse models are appropriate for some applications where VAR models are popularly adopted. Recently, Giannone et al. (2021) note the difficulty of identifying sparse predictive representations for several macroeconomic applications.
We illustrate the inadequacy of the sparsity assumption on a volatility panel dataset (see Section 5.3 for its description). Figure 1 (a) indicates that as the dimensionality increases, the leading eigenvalue of the spectral density matrix at frequency 0 (i.e. the long-run covariance) estimated from the data also increases linearly. This indicates the presence of strong serial and cross-sectional correlations that cannot be accommodated by sparse VAR models. Figure 1 (b) further confirms that we cannot infer a sparse pairwise relationship from such data. Once we estimate factors driving the strong correlations and adjust for their presence, however, there is evidence that the remaining dependence in the data can be modelled as being sparse, see Figure 1 (c). Together, the plots (d), (e), (c) and (f) display that the relationship between a pair of variables (after factor-adjustment) varies over time, particularly at the level of industrial sectors. Here, the intervals are chosen according to the data segmentation result reported in Section 5.3. This example highlights the importance of (i) accounting for the dominant serial and cross-sectional dependence prior to fitting a model under the sparsity assumption, and (ii) detecting possibly multiple change points in the underlying structure, when analysing time series datasets covering a long period.

Our contributions
Motivated by the aforementioned characteristics of high-dimensional time series data, factor-adjusted regression modelling has increasingly gained popularity (Fan et al., 2020(Fan et al., , 2021Krampe and Margaritella, 2021). The factor-adjusted VAR model proposed by Barigozzi et al. (2022) assumes that a handful of common factors capture strong serial and cross-sectional correlations, such that it is reasonable to assume a sparse VAR model on the remaining component to capture idiosyncratic, variable-specific dependence. We extend this framework by proposing a new, piecewise stationary factor-adjusted VAR model and develop FVARseg, an accompanying change point detection methodology. Below we summarise the methodological and theoretical contributions made in this paper.
Generality of the modelling framework. Under the proposed model, we decompose the data into two piecewise stationary latent processes: one is driven by factors and accounts for dominant serial and cross-sectional correlations, and the other models sparse pairwise dependence via a VAR model. We adopt the most general approach to factor modelling and allow both components to undergo changes which, in the case of the latter, are attributed to shifts in the VAR parameters. To the best of our knowledge, such a general model simultaneously permitting the presence of common factors and change points, has not been studied in the literature previously and provides a great deal of flexibility in high-dimensional time series modelling. Accordingly, we are not aware of any method that can comprehensively address the data segmentation problem considered in this paper.
Methodological novelty. The idea of scanning the data for changes over moving windows, has successfully been applied to a variety of data segmentation problems (Preuss et al., 2015;Yau and Zhao, 2016;Eichinger and Kirch, 2018;Chen et al., 2021). We propose FVARseg, a two-stage methodology that combines this idea with statistics carefully designed to have good detection power against different types of changes in the two latent components. In Stage 1 of FVARseg, motivated by the fact that dominant factor-driven correlations appear as leading eigenvalues in the frequency domain (see e.g. Figure 1 (a)), we propose a detector statistic that contrasts the local spectral density matrix estimators from neighbouring moving windows in operator norm, which is well-suited to detect changes in the factor-driven component.
In Stage 2 for detecting change points in the latent piecewise stationary VAR process, we deliberately avoid estimating the latent process which may incur large errors. Instead, we make use of (i) the Yule-Walker equation that relates autocovariances (ACV) and VAR parameters, and (ii) the availability of local ACV estimators of the latent VAR process after Stage 1. Combining these ingredients, we propose a novel detector statistic that enjoys methodological simplicity as well as statistical efficiency. Further, through sequential evaluation of the detector statistic, the second-stage procedure requires the estimation of local VAR parameters at selected locations only. Consequently it is highly competitive computationally when both the sample size and the dimensionality are large.
Theoretical consistency. FVARseg achieves consistency in estimating the total number and locations of the change points in both of the piecewise stationary factor-driven and VAR processes. Our theoretical analysis is conducted in a setting considerably more general than those commonly adopted in the literature, permitting dependence across stationary segments and heavy-tailedness of the data. We also derive the rate of localisation for each stage of FVARseg where we make explicit the influence of tail behaviour and the size of changes. In particular, under Gaussianity, the estimators from Stage 1 nearly matches the minimax optimal rate derived for the simpler, covariance change point detection problem, and Stage 2 also performs competitively with the state-of-the-art.
The rest of the paper is structured as follows. Section 2 introduces the piecewise stationary factor-adjusted VAR model. Section 3 describes the two stages of FVARseg, the proposed data segmentation methodology, and Section 4 establishes its theoretical consistency. Section 5 demonstrates the good performance of FVARseg on simulated datasets in comparison with competing methods, as well as applying it to a panel of volatility measures of US blue chip companies. R code implementing our method is available from https://github.com/haeran-cho/fvarseg.
Notation. Let I, O and 0 denote an identity matrix, a matrix of zeros and a vector of zeros whose dimensions depend on the context. For a random variable X and ν ≥ 1, denote X ν = (E|X| ν ) 1/ν . Given A = [a ii , 1 ≤ i ≤ m, 1 ≤ i ≤ n], we denote by A * its transposed complex conjugate. We define its element-wise ∞ , 1 and 2 -norms by |A| ∞ = max i,i |a ii |, |A| 1 = i,i |a ii | and |A| 2 = i,i |a ii | 2 , and its spectral and induced L 1norms by A and A 1 = max 1≤i ≤n m i=1 |a ii |, respectively. For positive definite A, we denote its minimum eigenvalue by A min . For two real numbers, a ∨ b = max(a, b) and a ∧ b = min(a, b). For two sequences {a n } and {b n }, we write a n b n if, for some constants C 1 , C 2 > 0, there exists N ∈ N such that C 1 ≤ a n b −1 n ≤ C 2 for all n ≥ N .
2 Piecewise stationary factor-adjusted VAR model

Background
A zero-mean, p-variate process ξ t follows a VAR(d) model if it satisfies where the parameter matrices A ∈ R p×p , 1 ≤ ≤ d, determine how future values of the series depend on their past. The p-variate random vector ε t = (ε 1t , . . . , ε pt ) has ε it which are independently and identically distributed (i.i.d.) for all i and t with E(ε it ) = 0 and Var(ε it ) = 1. Then, the positive definite matrix Γ ∈ R p×p is the covariance matrix of the innovations for the VAR process.
In the literature on factor modelling of high-dimensional time series, the factor-driven component exhibits strong cross-sectional and/or serial correlations by 'loading' a finitedimensional factors linearly. Among many models, the generalised dynamic factor model (GDFM, Forni et al., 2000 provides the most general approach (see Appendix C for further discussions) where the p-variate factor-driven component χ t admits For fixed q, the q-variate random vector u t = (u 1t , . . . , u qt ) contains the common factors which are shared across the variables and time, and u jt are assumed to be i.i.d. for all j and t with E(u jt ) = 0 and Var(u jt ) = 1. The matrix of square-summable filters B(L) = ∞ =0 B L with the lag-operator L and B ∈ R p×q , serves the role of loadings under (2). Barigozzi et al. (2022) propose a factor-adjusted VAR model for high-dimensional time series modelling and forecasting. Under the model, the observations are assumed to be decomposed as a sum of the two latent components ξ t and χ t in (1)-(2), whereby pervasive correlations in the data are accounted for by χ t and the remaining dependence is captured by ξ t . In the next section, we introduce its piecewise stationary extension where both the factor-driven and VAR processes are allowed to undergo structural changes.

Model
We observe a zero-mean, p-variate piecewise stationary process X t satisfying Here, θ χ,k , 1 ≤ k ≤ K χ , denote the change points in the piecewise stationary factor-driven component χ t such that at each θ χ,k , the filter of loadings B [k] (L) undergoes a change. We permit the factor number to vary over time as q k ≤ q, with the factor u t being a sub-vector of u t ∈ R q . Similarly, θ ξ,k , 1 ≤ k ≤ K ξ , denote the change points in the piecewise stationary VAR process ξ t at which the VAR parameters undergo shifts; we permit the VAR innovation covariance matrix to vary as Γ [k] but our interest lies in detecting changes in VAR parameters, and the VAR order may vary over time as By convention, we denote θ χ,0 = θ ξ,0 = 0 and θ χ,Kχ+1 = θ ξ,K ξ +1 = n. In line with the factor modelling literature, we assume that χ t and ξ t are uncorrelated through having E(u jt ε it ) = 0 for any i, j, t and t .
The model (3) does not require that the change points in χ t and ξ t are aligned, or that K χ = K ξ . Our goal is to estimate the total number and locations of the change points for both of the piecewise stationary latent processes. Importantly, we allow {ξ t , t ∈ Z}) to be dependent across k through sharing the innovations ε t (resp. u t ). This makes our model considerably more general than those found in the literature on (high-dimensional) data segmentation under VAR models (Wang et al., 2019;Safikhani and Shojaie, 2022;Bai et al., 2022) which assume independence across the segments.
Data segmentation under factor models has been considered by Barigozzi et al. (2018) and Li et al. (2022) where they adopt a static approach to factor modelling. On the other hand, the model in (3) adopts a piecewise stationary extension of GDFM which provides the most general approach to factor modelling, see Appendix C for further discussions.

Assumptions
In this section, we introduce assumptions that ensure the (asymptotic) identifiability of the two latent processes in (3) which are framed in terms of spectral properties, as well as controlling the degree of serial dependence in the data. We denote by Γ t at lag ∈ Z, and its spectral density matrix at frequency ω ∈ [−π, π] by Σ χ,j (ω), 1 ≤ j ≤ q k , denote the real, positive eigenvalues of Σ [k] χ (ω) ordered by decreasing size. We similarly define Γ Assumption 2.1. For each 0 ≤ k ≤ K χ , the following holds: There exist a positive integer p 0 ≥ 1, pairs of functions ω → α [k] j (ω) and ω → β [k] j (ω) for ω ∈ [−π, π] and 1 ≤ j ≤ q k , and r k,j ∈ (0, 1] satisfying r k,1 ≥ . . . ≥ r k,q k such that for all p ≥ p 0 , If r k,j = 1 for all 1 ≤ j ≤ q k as frequently assumed in the literature (Fan et al., 2013;Forni et al., 2015), we are in presence of q k factors which are equally pervasive for the whole cross-sections of χ [k] t . If r k,j < 1 for some j, we permit the presence of 'weak' factors. Since our primary interest lies in change point analysis, we later introduce a related but distinct condition on the size of change in χ t in Assumption 4.2.
Assumption 2.3. There exist constants Ξ > 0 and ς > 2 such that for all ≥ 0, Assumption 2.2 (i)-(ii) are standard conditions in the literature (Lütkepohl, 2005;Basu and Michailidis, 2015). Under condition (iii) and Assumption 2.3, we have serial dependence in X t decay at an algebraic rate according to the functional dependence measure of Zhang and Wu (2021). Assumption 2.2 (iii) allows for mild cross-sectional dependence in ξ [k] t while ensuring that µ [k] ξ,1 (ω) is uniformly bounded: Proposition 2.1. Under Assumption 2.2, uniformly over all ω ∈ [−π, π], there exists some M ξ > 0 depending only on M ε , Ξ and ς such that max 0≤k≤K ξ sup ω∈[−π,π] µ This, together with Assumption 2.2 (iv), establishes the boundedness of the eigenvalues of Σ ξ (ω), which is commonly assumed in the high-dimensional VAR literature for guaranteeing the consistency of Lasso estimators (Basu and Michailidis, 2015). Then from Assumption 2.1 and Proposition 2.1, the two latent components in (3) are asymptotically identifiable as p → ∞, thanks to the gap between µ [k] χ,q k (ω) diverging with p and µ [k] ξ,1 (ω) which is uniformly bounded. This also indicates that the proposed model (3) is able to capture the phenomenon observed in Figure 1 (a).

Methodology
In this section, we describe FVARseg, a two-stage data segmentation methodology which detects multiple change points in the factor-driven component χ t (Stage 1, Section 3.1) and the piecewise stationary VAR process ξ t (Stage 2, Section 3.2) under the model (3).

Change point detection
The spectral density matrix of χ t is given by Σ (e −ιω )) * for θ χ,k + 1 ≤ t ≤ χ χ,k+1 , i.e. it varies over time in a piecewise constant manner with change points at θ χ,k , 1 ≤ k ≤ K χ . By Weyl's inequality, Assumption 2.1 and Proposition 2.1 jointly indicate a gap in the eigenvalues of (time-varying) spectral density matrix of X t , i.e. those corresponding to the common component diverges with p while the remaining ones are bounded for all p. This suggests an approach that looks for changes in χ t from the behaviour of X t in the frequency domain which we further justify below.
Example 3.1. Suppose that χ t contains a single change point at t = θ χ,1 at which a new factor is introduced, i.e. χ . Let Σ x,t (ω) denote the time-varying spectral density of X t . Then, from the uncorrelatedness between χ t and ξ t and Proposition 2.1, That is, the change in the spectral density of χ t is detectable as a change in time-varying spectral density matrix of X t in operator norm, with the size of change diverging with p as (b(e −ιω )) * b(e −ιω ) does so under Assumption 2.1 on the strength of factors.
Thus motivated, we propose to detect and locate changes in χ t by scanning for any large change in the spectral density matrix of X t measured in operator norm. For this, we propose the following moving window-based methodology. Given a bandwidth G, we estimate the local spectral density matrix of X t by where K(·) denotes the Bartlett kernel, m = G β the kernel bandwidth with β ∈ (0, 1), and Then the following statistic serves as a good proxy of the difference in local spectral density matrices of χ t over I v (G) = {v − G + 1, . . . , v} and I v+G (G) = {v + 1, . . . , v + G}. To make it more precise, let Σ χ,v (ω, G) denote a weighted average , as a function of v, linearly increases and then decreases around the change points with a peak of size Σ (ω) formed at v = θ χ,k for all 1 ≤ k ≤ K χ , provided that the bandwidth G is not too large (in the sense of Assumption 4.2 (ii) below). The statistic T χ,v (ω, G) is designed to approximate T * χ,v (ω, G) when χ t is not directly observed, and thus is well-suited to detect and locate the change points therein. Algorithm 1 outlines the moving window-based procedure for multiple change point detection in χ t .
Step 1 of Algorithm 1, we evaluate T χ,v (ω, G) at the Fourier frequencies ω l = 2πl/(2m + 1) for 0 ≤ l ≤ m, and ω l = −ω −l for −m ≤ l ≤ −1, and take the maximum over the frequencies at each given location v, which then is compared against some threshold κ n,p . Collecting all time points at which max 0≤l≤m T χ,v (ω l , G) > κ n,p in I, there is some evidence for a change point around those time points belonging to I. However, care is needed to avoid detecting the same change point more than once, since T χ,v (ω l , G) is expected to take a large value over an interval containing θ χ,k for each k. Therefore in Step 2, we regard θ = arg max v∈I max l T χ,v (ω l , G) as a change point estimator if it is a local maximiser of T χ,v (ω( θ), G) within an interval of radius ηG centred at θ, with ω( θ) = arg max ω l : 0≤l≤m T θ (ω l , G) and some η ∈ (0, 1). Once θ is added to the set of final estimators, Step 3 ensures that we detect only a single change point for each θ χ,k by disregarding the interval of radius G centred at θ. By Step 4, we iteratively consider the maximisers of T χ,v (ω l , G) and determine whether they detect change points or not.

Post-segmentation factor adjustment
Following the detection of change points in χ t , we are able to estimate the segment-specific quantities related to χ [k] t . In view of the second-stage of FVARseg detecting change points in ξ t , we describe how to estimate Γ [k] χ ( ) with which we can estimate the ACV of ξ t . For each k = 0, . . . , K χ , we first estimate the spectral density of X t over the segment { θ χ,k +1, . . . , θ χ,k+1 } by Σ [k] x (ω) as in (4) using the sample ACV computed from the segment (we use the same kernel bandwidth m for simplicity). Then noting that the spectral density matrix of χ [k] t is of rank q k under (3), we estimate it from the eigendecomposition of Σ x (ω l ) by retaining only the q k largest eigenvalues, say µ [k] x,j (ω l ), and the associated eigenvectors e [k] x,j (ω l ), and the ACV of χ x,j (ω l ) e [k] x,j (ω l ) e [k] x,j (ω l ) * and Γ [k] The estimators in (7) require the factor number q k as an input. We refer to Hallin and Liška (2007) and Onatski (2009) for consistent estimators of q k that make use of the postulated eigengap in the spectral density matrix of X t .

Stage 2: Piecewise VAR process segmentation
Applying the existing VAR segmentation methods in our setting requires estimating the np elements of the latent piecewise stationary VAR process ξ t , which introduces additional errors and possibly results in the loss of statistical efficiency. In addition, as discussed in Appendix A.2, the existing methods tend to be computationally demanding, e.g. by evaluating the Lasso estimators O(n 2 ) times in a dynamic programming algorithm, or solving a large fused Lasso objective function of dimension np 2 d. Instead, since we can estimate the local AVC of ξ t from the post-segmentation factor-adjustment in Stage 1, we propose a methodology for segmenting the latent VAR component which avoids the need to estimate ξ t . As described below, the proposed method evaluates the local VAR parameters at carefully selected locations only, and thus is computationally efficient.
Specifically, our approach makes use of the Yule-Walker equation (Lütkepohl, 2005).
d ] ∈ R (pd)×p contain all VAR parameters in the kth segment. Then, it is related to the ACV matrices Γ being invertible due to Assumption 2.2 (iv). We propose to utilise this estimating equation in combination with the local ACV estimators of ξ t obtained as described below. For a given bandwidth G and the interval (7), with the weights given by the proportion of I v (G) covered by the kth segment (see (E.16) for the precise definition). Replacing Γ and similarly g v (G). Then, we propose to scan with some inspection parameter β ∈ R (pd)×p and a matrix norm ||| · |||. We motivate this statistic by considering T * if v is far from all the change points in ξ t , i.e. min k |v − θ ξ,k | ≥ G, while it is tent-shaped near the change points with a local maximum at v = θ ξ,k , provided that For the inspection parameter, we adopt an 1 -regularised Yule-Walker estimator of the VAR parameters first considered by Barigozzi et al. (2022) in stationary settings. At given v • ∈ {G, . . . , n}, we solve the constrained 1 -minimisation problem with a tuning parameter λ n,p > 0. The ∞ -constraint in (10) naturally leads to the choice ||| · ||| = | · | ∞ , resulting in the following detector statistic: For good detection power, the condition in (9) suggests using an estimator of β [k−1] or β [k] in place of β for detecting θ ξ,k . Therefore, we propose to evaluate computed at the three iterations of Steps 1-3 of Algorithm 2. At each iteration, we use β = β v• (G) estimated from each of the subsections the data highlighted in the x-axis (left to right); the corresponding estimators are plotted in the bottom panel and for comparison, we also plot the estimators obtained in the oracle setting where ξ t is observable (all plots have the identical z-axis range). The locations of v • , θ and θ in Algorithm 2, and θ ξ,k are denoted by the vertical long-dashed, dot-dashed, dotted and dashed lines, respectively, and the horizontal line represents π n,p chosen as described in Section 5.1. (10) with v • = G and scan the data using When T ξ,v ( β, G) exceeds some threshold, say π n,p , at v =θ for the first time, it signifies that a change has occurred in the neighbourhood. Reducing the search for a change point to {θ, . . . ,θ + G}, we identify a change point estimator as the local maximiser θ ξ,1 = arg maxθ ≤v≤θ+G T ξ,v ( β, G). Then updating β with β v• (G) obtained at v • = θ ξ,1 + (η + 1)G for some η ∈ (0, 1] (i.e. only using an interval of length G located strictly to the right of θ ξ,1 for its estimation), we continue screening T ξ,v ( β, G), v ≥ v • , until it next exceeds π n,p . These steps of screening T ξ,v ( β, G) and updating β are repeated iteratively until the end of the data sequence is reached.
Algorithm 2 outlines the steps of the Stage 2 methodology, and Figure 2 provides an illustration. In particular, it shows that although ξ t is not observed, at each iteration, β v• (G) does as well as its oracle counterpart (obtained as in (10) with the sample ACV of ξ t replacing Γ ξ,v ( , G)). Computationally, this strategy benefits from that the costly Algorithm 2: Stage 2 of FVARseg Input: Data {X t } n t=1 , change point estimators from the common component Θ χ , λ n,p for (10), bandwidth G, η ∈ (0, 1], threshold π n,p Step 0: Set Θ ξ ← ∅ and v • ← G.
Output: Θ ξ solution to the 1 -minimisation problem in (10) is required (at most) K ξ + 1 times with an appropriately chosen threshold π n,p (see Theorem 4.3 below). We further demonstrate numerically the computational competitiveness of Stage 2 as a standalone method for VAR time series segmentation in Section 5.2, and provide an in-depth comparative study with the existing methods in Appendix A.2.

Consistency of Stage 1 of FVARseg
We carry out our theoretical investigation under two different regimes with respect to the tail behaviour of u t and ξ t ; in particular, the weaker condition in Assumption 4.1 (i) permits heavy-tailed innovations, while the existing literature on (piecewise stationary) VAR modelling in high dimensions, commonly adopts the Gaussianity as in (ii).
Assumption 4.1. We assume either of the following conditions.
In establishing the consistency of Stage 1, which detects changes that are attributed to changes in the factor structure, we opt to measure the size of changes in χ t using the difference in spectral density matrices of χ t from neighbouring segments. As ∆ χ,k (ω) is Hermitian, we can always find the jth largest (in modulus), real-valued eigenvalue of ∆ χ,k (ω) which we denote by µ j (∆ χ,k (ω)), with µ 1 (∆ χ,k (ω)) = ∆ χ,k (ω) . Recall that m = G β for some β ∈ (0, 1), denotes the bandwidth used in local spectral density estimation, see (4).
Assumption 4.2. (i) For each 1 ≤ k ≤ K χ , the following holds: There exist a positive integer p 0 ≥ 1 and pairs of functions ω → a [k] j (ω) and ω → b [k] j (ω) for ω ∈ [−π, π] and j = 1, 2, and r k,1 ∈ (0, 1] and r k,2 ∈ [0, 1] satisfying r k,1 ≥ r k,2 , such that Assumption 4.2 specifies the detection lower bound which is determined by min k ∆ χ,k and min k (θ χ,k+1 −θ χ,k ) (through G), for all K χ change points χ t to be detectable by Stage 1. Condition (i) requires µ 1 (∆ χ,k (ω)) to be distinct from the rest. In fact, the remaining µ j (∆ χ,k (ω)), j ≥ 2, are allowed to be exactly zero, which is the case in Example 3.1; here, we have ∆ χ,1 = max ω (2π) −1 (b(e −ιω )) * b(e −ιω ) where b(z) is a p-variate vector of factor loading filters. The rate p(ψ n ∨ m −1 ) represents the bias-variance trade-off when estimating the local spectral density matrix of χ t by Σ x,v (G, ω) (see Proposition E.6). It is possible to find the rate of kernel bandwidth m that minimises this rate depending on the tail behaviour of X it (e.g. m (G/ log(n)) 1/3 under Gaussianity), but we choose to explicitly highlight the role of this tuning parameter on our results.
n,p reflects the difficulty associated with estimating the individual change point θ χ,k manifested by (p −1 ∆ χ,k ) −2 . In the Gaussian case (Assumption 4.1 (ii)), the localisation rate ρ [k] n,p is always sharper than G due to Assumption 4.2 (i). Considering the problem of covariance change point detection in independent, sub-Gaussian random vectors in high dimensions, Wang et al. (2021) derive the minimax lower bound on the localisation rate in their Lemma 3.2, and ρ [k] n,p matches this rate up to m log(n); here, the dependence on the kernel bandwidth m is attributed to that we consider a time series segmentation problem, i.e. a change may occur in the ACV of χ t at lags other than zero. If heavier tails are permitted (Assumption 4.1 (i)), ρ [k] n,p can be tighter than 0 G e.g. when ∆ χ,k p, K χ is fixed and m G β for some β ∈ (0, 1 − 4/ν).
(ii) Empirically, replacing θ with θ = arg max v∈I avg l T χ,v (ω l , G) in Step 2 of Algorithm 1 returns a more stable location estimator, where avg l denotes the average operator over l = 0, . . . , m. We can derive the localisation rate for θ similar to that reported in Theorem 4.1 (b) with ∆ χ,k = π −1 π 0 ∆ χ,k (ω) dω in place of ∆ χ,k , and our numerical results in Section 5.2 are based on this estimator.
Next, we establish the consistency of Γ t under the following assumption on the strength of factors.
Theorem 4.2. Suppose that Assumption 4.3 holds in addition to the assumptions made in Theorem 4.1, and define ρ n,p = max 1≤k≤Kχ min( 0 G, ρ [k] n,p ). Also let Then on M χ n,p defined in Theorem 4.1, for some finite integer d ∈ N, we have It is possible to work under the weaker Assumption 2.1 and trace the effect of weak factors or bound estimation errors measured in different norms, as in Barigozzi et al. (2022) deriving such results in the stationary setting. We choose to work under the stronger Assumption 4.3 as it simplifies the presentation of Theorem 4.2 which plays an important role in the investigation into Stage 2 of FVARseg given in Section 4.2 below.

Consistency of Stage 2 of FVARseg
Suppose that the tuning parameter for the 1 -regularised Yule-Walker estimation problem in (10), is set with some constant M > 0 and ϑ n,p defined in (13), as This choice reflects the error in Γ ξ,v ( , G) estimating the local ACV of ξ t over all v and .
The following assumption imposes conditions on the size of the changes in VAR parameters and the minimum spacing between the change points.
Remark 4.2. We choose to measure the size of change using | 2 scaled by the global sparsity (given by the union of the supports of all , is used to measure the size of change where this global sparsity may be much greater than that of ∆ ξ,k when K ξ is large. When d = 1, if the lag-zero ACV of ξ t remains constant before and after the change point θ ξ,k , we have G [k] (G [k−1] ) −1 1 = 1 and Assumption 4.4 (i) becomes min k |∆ ξ,k | ∞ /λ n,p → 0. More generally, bounding 1 implicitly assumes (approximate) sparsity on the second-order structure of ξ t . Han et al. (2015) show that if all G [k] are strictly diagonally dominant (see Definition 6.1.9 of Horn and Johnson, 1985) with min 1≤i≤p (|γ Theorem 4.3. Suppose that Assumption 4.4 holds in addition to the assumptions made in Theorem 4.2. With λ n,p chosen as in (14), we set π n,p to satisfy 2λ n,p < π n,p < 1 2 min Then, there exists a set M ξ n,p with P(M ξ n,p ) → 1 as n, p → ∞, such that the following holds for Θ ξ = { θ ξ,k , 1 ≤ k ≤ K ξ : θ ξ,1 < . . . < θ ξ, K ξ } returned by Stage 2 of FVARseg, on M ξ n,p for large enough n: under Assumption 4.1 (ii).
Theorem 4.3 (a) establishes that Stage 2 of FVARseg consistently detects all K ξ change points within the distance of 0 G where 0 can be made arbitrarily small as n, p → ∞ under Assumption 4.4 (i). Theorem 4.3 (b) shows that a further refined localisation rate can be derived for θ ξ,k when it is sufficiently distanced away from the change points in the factor-driven component. This is due to the sequential nature of FVARseg where the second-stage depends on the output from the first. If, say, θ ξ,k lies close to θ χ,k , a change point in χ t , the error from estimating the local ACV of ξ t due to the bias in θ χ,k , prevents applying the arguments involved in the refinement to such θ ξ,k . The refined rate [k] n,p is always tighter than G under Gaussianity.
It is of independent interest to consider the cases where the factor-driven component χ t is stationary (i.e. K χ = 0) or where we directly observe the piecewise stationary VAR process (i.e. X t = ξ t ). Consistency of the Stage 2 methodology of FVARseg readily extends to such settings and the improved localisation rates in Theorem 4.3 (b) apply to all the estimators. Also, further improvement is attained in the heavy-tailed situations (Assumption 4.1 (i)) in the case where we directly observe ξ t . For the full statement of the theoretical results, we refer to Corollary A.1 in Appendix A where we also provide a detailed comparison between Stage 2 of FVARseg and existing VAR segmentation methods (that do not take into the possible presence of factors), both theoretically and numerically.

Numerical considerations
Multiscale extension. The bandwidth G is required to be large enough to provide a good local estimators of spectral density of χ t (Stage 1) and VAR parameters (Stage 2). However, if G is too large, we may have windows that contain two or more changes when scanning the data for change points, which violates Assumptions 4.2 (ii) and 4.4 (ii). Cho and Kirch (2022) note the lack of adaptivity of a single-bandwidth moving window procedure in the presence of multiscale change points (a mixture of large changes over short intervals and smaller changes over long intervals), and advocates the use of multiple bandwidths. Accordingly we also propose to apply FVARseg with a range of bandwidths and prune down the outputs using a 'bottom-up' method of Messer et al. (2014); Meier et al. (2021). Let Θ(G) denote the output from Stage 1 or 2 with a bandwidth G. Given a set of bandwidths G = {G h , 1 ≤ h ≤ H : G 1 < . . . < G H }, we accept all estimators from the finest G 1 to the set of final estimators Θ and sequentially for h ≥ 2, accept θ ∈ Θ(G h ) iff minθ ∈ Θ | θ −θ| ≥ G/2. In simulation studies, we use G χ = {[n/10], [n/8], [n/6], [n/4]} for Stage 1, and G ξ generated as an equispaced sequence between [2.5p] and [n/4] of length 4 for Stage 2. The choice of G ξ is motivated by the simulation results of Barigozzi et al. (2022) under the stationarity, where the 1 -regularised estimator in (10) was observed to performs well when the sample size exceeds 2p.
Speeding up Stage 1. The computational bottleneck of FVARseg is the computation of T χ,v (ω l , G) in Stage 1, which involves singular value decomposition (SVD) of a p × pmatrix at multiple frequencies and over time. We propose to evaluate T χ,v (ω l , G) on a grid v ∈ {G + ab n : 0 ≤ a ≤ (n − 2G)/b n } with b n = 2 log(n) . This may incur additional bias of at most b n /2 ≤ log(n) in change point location estimation which is asymptotically negligible in view of Theorem 4.1, but reduce the computational load by the factor of b n .
Selection of thresholds. The theoretically permitted ranges of κ n,p and π n,p (see Theorems 4.1 and 4.3) depend on constants which are not accessible or difficult to estimate in practice. This is an issue commonly encountered by data segmentation methods involving localised testing, and often a reasonable solution is found by large-scale simulations, an approach we also take. We simulate B = 100 time series following (3) with K χ = K ξ = 0 using the models considered in our simulation studies (Section 5.2), and record the maximum of where the scaling terms are obtained from the first G observations only, with the aim of having the resultant statistics less dependent on the data generating process. Repeating the same procedure for varying (n, p, G), we fit a linear model to the 100(1−τ )th percentile of log(max v T • χ,v (G)) with log log(n) and log(G) as regressors (R 2 adj = 0.9651), and use the fitted model to derive a threshold for given n and G that is then applied to the similarly scaled T χ,v (ω l , G). Analogously, we regress the 100(1 − τ )th percentile of log(max v T • ξ,v (G)) onto log log(n), log log(p) and log(G) (R 2 adj = 0.985), and find a threshold applied to the scaled T ξ,v ( β, G) given n, p and G from the fitted model. When Stage 2 is used as a standalone method for segmenting observed VAR processes, a smaller threshold is recommended which is in line with Corollary A.1, and we find that π n,p = 1 works well with the proposed scaling.
Selection of other tuning parameters. There are data-adaptive methods for selecting the kernel window size m in (4) that work well in finite dimensions (Politis, 2003). Since our aim lies in change point analysis rather than in accurate estimation of the spectral density matrix, we propose to set it simply at m = max(1, G 1/3 ) for given G. The results are not highly sensitive to the choice of η in Stage 1 and use η = 0.5 throughout. In Stage 2, we find that not trimming off the data when estimating the VAR parameters by setting η = 0, does not hurt the numerical performance. In factor-adjustment, we select the segment-specific factor number q k using the information criterion (IC)-based approach of Hallin and Liška (2007). Krampe and Margaritella (2021) propose to jointly select the (static) factor number and the VAR order under their factor-adjusted VAR model using an IC but generally, the validity of IC is not well-understood for the task of VAR order selection in high dimensions. In our simulations, following the practice in the VAR segmentation literature, we regard d as known while in analysing the panel of daily volatilities (Section 5.3), we use d = 5 which has the interpretation of the number of trading days per week. Finally, we select λ n,p in (10) via cross-validation as in Barigozzi et al. (2022).

Simulation studies
In the simulations, we consider the cases when the factor-driven component is present (χ t = 0) and when it is not (χ t = 0). For the former, we consider two models for generating χ t with q = 2. In the first model, referred to as (C1), χ t admits a static factor model representation while in the second model (C2), it does not; empirically, the task of factor structure estimation is observed to be more challenging under (C2) (Forni et al., 2017;Barigozzi et al., 2022). Full descriptions of (C1) and (C2) can be found in Appendix B.1. For generating the piecewise stationary VAR(d) process ξ t , we consider  Table 1 for an overview of the 24 data generating processes used in our simulation studies which also contains information about the sets of change points Θ χ and Θ ξ ; under each setting, we generate 100 realisations. Below we provide a summary of the findings from the simulation studies, and Tables B.1-B.2 reporting the results can be found in Appendix B.2. To the best of our knowledge, there does not exist a methodology that comprehensively addresses the change point problem under the model (3). Therefore under (M1)-(M2), we compare the Stage 1 of FVARseg with a method proposed in Barigozzi et al. (2018), referred to as BCF hereafter, on their performance at detecting changes in χ t . While BCF has a step for detecting change points in the remainder component, it does so nonparametically unlike the Stage 2 of FVARseg, which may lead to unfair comparison. Hence we separately consider (M3) with X t = ξ t where we compare the Stage 2 method with VARDetect, a block-wise variant of the method proposed in Safikhani and Shojaie (2022), see also  and Bai et al. (2021).
Results under (M1)-(M2). Overall, FVARseg achieves good accuracy in estimating the total number and locations of the change points for both χ t and ξ t across different data generating processes. Under (M1) adopting the static factor model for generating χ t , FVARseg shows similar performance as BCF in detecting Θ χ when the dimension is small (p = 50), but the latter tends to over-estimate the number of change points as p increases. Also, FVARseg outperforms the binary segmentation-based BCF in change point localisation. BCF requires as an input the upper bound on the number of global factors, say q , that includes the ones attributed to the change points, and its performance is sensitive to its choice. In (M1), we have q ≤ 3q(K χ + 1) (which is supplied to BCF) while in (M2), χ

[k]
t does not admit a static factor representation and accordingly such q does not exist (we set q = 2q for BCF). Accordingly, BCF tends to under-estimate the number of change points under (M2). Generally, the task of detecting change points in ξ t is aggravated by the presence of change points in χ t due to the sequential nature of FVARseg, and the Stage 2 performs better when K χ = 0 both in terms of detection and localisation accuracy, which agrees with the observations made in Corollary A.1 (a). Between (M1) and (M2), the latter poses a more challenging setting for the Stage 2 methodology. This may be attributed to (i) the difficulty posed by the data generating scenario (C2), which is observed to make the estimation tasks related to the latent VAR process more difficult (Barigozzi et al., 2022), and (ii) that Θ χ = Θ ξ where the estimation bias from Stage 1 has a worse effect on the performance of Stage 2 compared to when Θ χ and Θ ξ do not overlap, see the discussion below Theorem 4.3.
Results under (M3). Table B.2 shows that the Stage 2 of FVARseg outperforms VARDetect in all criteria considered including the execution time, particularly as p increases. VARDetect struggles to detect any change point when the change is weak (recall that β = 0.6 is used when d = 1 which makes the size of change at θ ξ,2 small) or when d = 2. FVARseg is faster than VARDetect in all situations, sometimes up to 10 times e.g. when d = 2 and there is no change point in the data. Compared to the results obtained under (M1)-(M2), the localisation performance of the Stage 2 method improves in the absence of the factor-driven component, even though the size of changes under (M3) tends to be smaller; this confirms the theoretical findings reported in Corollary A.1 (b) in Appendix A. Although not reported here, when the full FVARseg methodology is applied to the data generated under (M3), the Stage 1 method does not detect any spurious change point estimators as desired.

Application: US blue chip data
We consider daily stock prices from p = 72 US blue chip companies across industry sectors between January 3, 2000 and February 16, 2022 (n = 5568 days), retrieved from the Wharton Research Data Services; the list of companies and their corresponding sectors can be found in Appendix D. Following Diebold and Yılmaz (2014), we measure the volatility using (resp. p low it ) denotes the maximum (resp. minimum) log-price of stock i on day t, and set X it = log(σ 2 it ). We apply FVARseg to detect change points in the panel of volatility measures {X it , 1 ≤ i ≤ p; 1 ≤ t ≤ n}. With n 0 = 252 denoting the number of trading days per year, we apply Stage 1 with bandwidths chosen as an equispaced sequence between [n 0 /4] and 2n 0 of length 4, implicitly setting the minimum distance between two neighbouring change points to be three months. Based on the empirical sample size requirement for VAR parameter estimation (see Section 5.1), we apply Stage 2 with bandwidths chosen as an equispaced sequence between 2.5p and 2n 0 of length 4. The VAR order is set at d = 5 which corresponds to the number of trading days in each week, and the rest of the tuning parameters are selected as in Section 5.1. Table 2 reports the segmentation results.   To further validate the change point estimators obtained by FVARseg, we perform a forecasting exercise. Two approaches, referred to as (F1) and (F2) below, are adopted to build forecasting models where the difference lies in how a sub-sample of {X u , u ≤ t−1}, is chosen to forecast X t , t ∈ T ⊂ {1, . . . , n}. Simply put, (F1) uses the observations belonging to the same segment as t only, for constructing the forecast of χ t (resp. ξ t ) according to the segmentation defined by Θ χ (resp. Θ ξ ), while (F2) ignores the presence of the most recent change point estimator. We expect (F1) to give more accurate predictions if the data undergoes structural changes at the detected change points. On the other hand, if some of the change point estimators are spurious, (F2) is expected to produce better forecasts since it makes use of more observations. More precisely, defining L χ (v) = max{0 ≤ k ≤ K χ : θ χ,k +1 ≤ v} denoting the index of θ χ,k nearest to and strictly left of v and similarly L ξ (v), we select T such that each t ∈ T fulfils min( L χ (t), L ξ (t)) ≥ 2 and min(t − θ χ, Lχ(t) , t − θ ξ, L ξ (t) ) ≥ n 0 ; this ensures that there are at least n 0 observations to produce a forecast. We have |T | = 1600. For such t ∈ T , we obtain X t (N) = χ t (N 1 ) + ξ t (N 2 ) for some N = (N 1 , N 2 ), where χ t (N 1 ) denotes an estimator of the best linear predictor of χ t given X t− , 1 ≤ ≤ N 1 , and ξ t (N 2 ) is defined analogously. We select (N 1 , N 2 ) using the following two approaches. Barigozzi et al. (2022) propose two methods for estimating the best linear predictors of χ t and ξ t under a stationary factor-adjusted VAR model, one based on a more restrictive assumption on the factor structure ('restricted') than the other ('unrestricted'); we refer to the paper for their detailed descriptions. Both estimators are combined with the two approaches (F1) and (F2). Table 3 reports the summary of the forecasting errors measured as FE avg t = |X t − X t (N)| 2 2 /|X t | 2 2 and FE max t = |X t − X t (N)| ∞ /|X t | ∞ , obtained from combining different best linear predictors with (F1)-(F2). According to all evaluation criteria, (F1) produces forecasts that are more accurate than (F2) regardless of the forecasting methods, which supports the validity of the change point estimators returned by FVARseg.

A.1 Extension of Theorem 4.3
We consider the performance of the Stage 2 of FVARseg when applied to some special cases under the model (3) where (a) χ t is stationary (i.e. Θ χ = ∅) and (b) we directly observe X t = ξ t (i.e. χ t = 0). Theorem 4.1 indicates that in both cases, the Stage 1 of FVARseg returns Θ χ = ∅. The results reported in Theorem 4.3 readily extend to such settings.
Corollary A.1. Suppose that the assumptions of Theorem 4.3 hold, including Assumption 4.4 (i) with λ n,p specified below. Then, with [k] n,p defined as in Theorem 4.3, i.e. [k] there exist a set M ξ n,p with P(M ξ n,p ) → 1 as n, p → ∞ and constants 0 , c 0 > 0 such that on M ξ n,p , we have n,p for all 1 ≤ k ≤ K ξ for n large enough, in the following situations.
(a) There is no change point in the factor-driven component, i.e. K χ = 0, and we set (b) We directly observe the piecewise stationary VAR process, i.e. X t = ξ t for all t, and we set λ n,p = M (max 0≤k≤K ξ β [k] 1 + 1)θ n,p with ( When compared to the methods dedicated to the setting corresponding to Corollary A.1 (b), our Stage 2 methodology achieves comparative theoretical performance in terms of the detection lower bound imposed on the size of changes for their detection, and the rate of localisation achieved. We provide a comprehensive comparison of the Stage 2 methodology with the existing VAR segmentation methods in the next section, both on their theoretical and computational properties.

A.2 Comparison with the existing VAR segmentation methods
There are a few methods proposed for time series segmentation under piecewise stationary, Gaussian VAR models, a setting that corresponds to Corollary A.1 (b) under Gaussianity. In this setting, we compare the Stage 2 of FVARseg with those proposed by Wang et al. (2019) and Safikhani and Shojaie (2022).
∈ S for all k. LP(a) denotes the complexity of solving a linear program with a variables and Lasso(a, b) that of solving a Lasso problem with sample size a and dimensionality b.

Methods
Separation Localisation Complexity Table A.1 summarises the comparative study in terms of their theoretical and computational properties. Denoting by∆ k the size of change between the (k − 1)th and the kth segments (measured differently for different methods), the separation rate refers to some ν n,p → ∞ such that if ν −1 n,p min k∆ 2 k · min k (θ ξ,k+1 − θ ξ,k ) → ∞, the corresponding method correctly detects all K ξ change points; for Stage 2, we set∆ k = |∆ ξ,k | ∞ and for the others, Table A.1 for the definition of s • ). The localisation rate refers to some n,p → ∞ satisfying max 1≤k≤K ξ w 2 k |θ k − θ ξ,k | = O P ( n,p ) for the estimatorsθ k returned by respective methods. For FVARseg, the weights w k reflect the difficulty associated with locating individual change points, i.e. w k =∆ k , while for Wang et al. (2019) and Safikhani and Shojaie (2022), the weights are global with w k = min k∆k and w k = 1, respectively. Safikhani and Shojaie (2022) further assume that min k |β [k] − β [k−1] | 2 is bounded away from zero. We suppose that max k β [k] 1 = O(1), a sufficient condition for the boundedness of Σ [k] ξ (ω) for each segment-specific VAR process (Basu and Michailidis, 2015, Proposition 2.2), which is required by all the methods in consideration for their theoretical consistency.
Immediate comparison of the theoretical results is difficult due to different definitions of∆ k : Observe that from Assumption 2.2 (iv), where | · | 0 denotes the element-wise 0 -norm. Noting that  Safikhani and Shojaie (2022) is generally worse than the other two both in terms of separation and localisation rates. The 1 -regularised Yule-Walker estimation problem in (10) can be solved in parallel and further, it needs to be performed only K ξ + 1 times with large probability, which makes the Stage 2 methodology more attractive. By comparison, the dynamic programming methodology of Wang et al. (2019) requires the Lasso estimation to be performed O(n 2 ) times, and the multi-stage procedure of Safikhani and Shojaie (2022) solves a fused Lasso problem of dimension np 2 d to obtain pre-estimators of the change points, and then exhaustively searches for the final set of estimators which can be NP-hard in the worst case. In Section 5.2, we compare the Stage 2 methodology with a blockwise modification of Safikhani and Shojaie (2022) that is implemented in the R package VARDetect (Bai et al., 2021).
Finally, we note that there are methods developed under piecewise stationary extensions of the low-rank plus sparse VAR(1) model proposed in Basu et al. (2019), see Bai et al. (2020Bai et al. ( , 2022. While they additionally permit a low rank structure in the parameter matrices, the spectrum of X t is assumed to be uniformly bounded which rules out pervasive (serial) correlations in the data and thus is distinguished from the piecewise stationary factor-adjusted VAR model considered in this paper.

B.1 Data generating processes for the common component
We provide full details on how the common component is generated for numerical experiments reported in Section 5.2.
t admits a static factor model representation, as where u jt ∼ iid N (0, σ 2 j ) with (σ 1 , σ 2 ) = (1, 0.5), and the MA coefficients are generated as (B t does not admit a static factor model representation, as where u jt ∼ iid N (0, 1) and the coefficients a ij are drawn uniformly as a ij ∼ iid U[− χ .

B.2 Complete simulation results
Tables B.1 and B.2 report the complete results obtained for the simulation studies described in Section 5.2. In particular, Table B.1 compares the performance of FVARseg against BCF (Barigozzi et al., 2018) on datasets generated as in (M1)-(M2) of Table 1, and Table B.2 compares the Stage 2 methodology of FVARseg (i.e. Algorithm 2 applied with Θ χ = ∅ and Γ ξ,v ( , G) = Γ x,v ( , G)), against VARDetect (Safikhani and Shojaie, 2022;Bai et al., 2021) on datasets generated under (M3) in Table 1. All tuning parameters are selected as described in Section 5.1. Denoting by Θ and Θ the sets of estimated and true change points, respectively, we report the distributions of K − K (with K = | Θ| and K = |Θ|) and the (scaled) Hausdorff distance between Θ and Θ, averaged over 100 realisations, as well as the average computation time (in seconds) in Table B.2. Table B.1: (M1)-(M2): Distributions of K χ − K χ and K ξ − K ξ and the average Hausdorff distance d H ( Θ χ , Θ χ ) and d H ( Θ ξ , Θ ξ ) returned by FVARseg and BCF (Barigozzi et al., 2018), over 100 realisations. We have K ξ = 2 under both (M1) and (M2).   (Bai et al., 2021), over 100 realisations. We also report the average computation time (in seconds) from 16 cores of an Intel Xeon Gold 6248R 3.00 GHz CPU with 16 GB of RAM on Linux (for d = 1) and 10 cores of Apple M1 Max with 64 GB of RAM on mac OS (for d = 2). C Generalised dynamic factor model C.1 GDFM as a representation Forni and Lippi (2001) show that the necessary and sufficient condition for any p-dimensional, weakly stationary time series to admit the generalised dynamic factor model (GDFM) representation, is to have a finite number of the eigenvalues of its spectral density matrix diverge with p (as in Assumption 2.1) while the remaining ones are bounded for all p. In other words, GDFM itself (without the VAR model imposed on ξ t as in this paper) can be regarded as a representation of high-dimensional time series rather than a model. Overall, GDFM provides the most general framework for high-dimensional time series factor modelling and it encompasses other factor models found in the literature such as static factor models (Forni et al., 2009). Static factor models are popularly adopted in both stationary (Stock and Watson, 2002;Bai, 2003;Fan et al., 2013;Barigozzi and Cho, 2020) and piecewise stationary (Barigozzi et al., 2018;Li et al., 2022) time series modelling in high dimensions. Under stationary factor models, the factor-driven component permits a representation χ it = λ i f t with some finite-dimensional vector processes f t ∈ R r as the common factors; here, 'static' refers to that χ it loads f t contemporaneously and does not preclude serial dependence therein. The model in (2) includes such a static factor model by representing Forni et al. (2009)). On the other hand, some models that have a finite number of factors under (2) cannot be represented with f t of finite dimension, the simplest example being the case where χ it = a i (1 − b i L) −1 u t for some b i ∈ (−1, 1) (Forni et al., 2015); see also (C2) in Section 5.2. Hallin et al. (2018) observe that principal component analysis (PCA), typically accompanying static factor models as an estimation tool, does not enjoy the optimality property that guarantees their success in the i.i.d. case in the presence of serial correlations, unlike the dynamic PCA adopted for estimation under GDFM (see Section 3.1.2).

C.2 VAR representation of GDFM
For notational simplicity, let q k = q for all 0 ≤ k ≤ K χ . Suppose that each (i, j)th element of the filter ,ij L , is a ratio of finite-order polynomials in L such that for some finite s 1 , s 2 ∈ N, ,ij L , l = 1, 2, for all 1 ≤ j ≤ q k and 0 ≤ k ≤ K χ . Furthermore, assume the followings.
(a) There exists M χ > 0 such that Under such assumptions, Section 4 of Forni et al. (2015) establishes that for generic values of the parameters B [k,1] ,ij and B [k,2] ,ij (outside a countable union of nowhere dense subsets), χ t admits a block-wise singular VAR representation 2] . . .
∈ R p×q is of rank q k ; for convenience, we assume that p = N (q + 1) for some N ∈ N. Here, each χ χ, L with its degree s ≤ qs 1 + q 2 s 2 , and det(A [k,h] χ (z)) = 0 for all |z| ≤ 1.
The representation (C.1) gives the piecewise stationary factor-adjusted VAR model in (3) the interpretation of decomposing high-dimensional time series into two latent VAR processes with time-varying parameter matrices, one of low rank (singular) accounting for dominant dependence and the other modelling individual interdependence between the variables unaccounted for by the former. Table D.1 provides the list of the 72 companies included in the application presented in Section 5.3 along with their tickers and industry classifications . (L)u t even when q k < q.
Proof. By Minkowski inequality, for all t. Due to independence of u jt , Assumption 2.3 and Lemma D.3 of Zhang and Wu (2021), there exists C ν > 0 that depends only on ν such that Similarly, from Assumption 2.2 and independence of ε it , we have and setting α ≤ ς − 1, Lemma E.2. There exist some constants C Ξ,ς , C Ξ,ς,ε > 0 which depend only on Ξ, ς and M ε defined in Assumptions 2.2 and 2.3, such that for all h ≥ 0, Proof. Suppose that L χ (t − ) = k and L χ (t) = l. Then, for any h ≥ 0, we have uniformly in 1 ≤ i, i ≤ p and t for some C Ξ,ς > 0 depending only on Ξ and ς, thanks to Assumption 2.3. Similarly, assuming that L ξ (t − ) = k and L ξ (t) = l, we have uniformly in 1 ≤ i, i ≤ p and t for some C Ξ,ς,ε > 0 depending only on Ξ, ς and M ε , from Assumption 2.2 (ii) and (iii).
The following lemma is a direct consequence of Lemma E.2.
such v cannot be a local maximiser of T χ,v (ω(v), G) within its ηG-radius. This, combined with how I is updated at each iteration, makes sure that only a single estimator is added to Θ χ for each change point.
Proof of Theorem 4.1 (b). WLOG, we consider the case when θ k ≤ θ k ; the following arguments apply analogously to the case when θ k > θ k . We prove by contradiction that if n,p , we have T χ, θ k (ω( θ k ), G) < T χ,θ k (ω( θ k ), G) and thus θ k cannot be the local maximiser of T χ,v (ω( θ k ), G) within its ηG-environment as required. From , G) and by (E.2) and (E.3), we have for an arbitrarily small constant > 0. Noting that Σ x,v (ω, G) and Σ χ,v (ω, G) are Hermitian (and thus diagonalisable with real diagonal entries), we write for some g k , g k ∈ C p satisfying g k = g k = 1, for all 1 ≤ k ≤ K χ . By Assumption 4.2 (ii), we have Then on E (1) n,p , there exists s k with |s k | = 1 and a constant c 1 > 0 such that where the first inequality follows from Corollary 1 of Yu et al. (2015), and the second one First, note that by (E.4), and as applying the same arguments as those adopted in bounding I and II above, we have max(R k1 , R k2 ) ≤ |F k |. Now we turn our attention to R k4 . Let γ x,ii (t, ) = γ By Lemma E.2, there exists constant C Ξ,ς,ε , C Ξ,ς,ε > 0 that do not depend on i, i such that Similarly, noting that there are at most a single change point within any (θ k − θ k )-interval, we have n,p . Collecting the bounds on IV , V and V I, we have for some constant C > 0, under Assumption 4.2 (i). Also, we observe that by Proposition 2.1, Turning our attention to R k3 , note that where the definitions of Q (r) k,ii (ω, h, H), r = 1, 2, can be found in (E.11). Then by Lemma E.8 and Chebyshev's inequality, there exists some constant c 1 > 0 such that P(Ē for some 1 ≤ δ ≤ G, where w k = (p −1 ∆ χ,k ) −1 and I k is defined in the lemma. Setting n,p (which itself does not depend on k), we have onĒ for c 0 large enough which, combined with the bounds on R kl , l = 1, 2, 4, 5, contradicts the first inequality in (E.7). As these statements are deterministic on E (1) n,p ∩Ē (1) n,p , the above arguments apply to all 1 ≤ k ≤ K χ which concludes the proof.

E.3.1 Supporting results
In what follows, we operate under the assumptions made in Theorem 4.1. (11).
Proof. Noting that we first address the first term in the RHS of (E.9). In Lemma E.1, for ς > 2 (as assumed in Assumptions 2.2 (iii) and 2.3), we can always set α = ς − 1 > 1/2 − 2/ν. Then, from the finiteness of Φ ν,α shown therein and by Theorems 4.1 and 4.2 of Zhang and Wu (2021), there exist universal constants C 1 , C 2 , C 3 > 0 and constants C α , C ν,α > 0 that depend only on their subscripts, such that for any z > 0, Noting that for any positive random variable Y , we have n for some constant C > 0 independent of i, i , thanks to Lemma E.1.
Turning our attention to the second term in the RHS of (E.9), let Γ Then by Lemma E.2, for all , 1 ≤ i, i ≤ p and G ≤ v ≤ n − G, From the bounds on I and II (which hold uniformly over 1 ≤ i, i ≤ p and G ≤ v ≤ n) and that ς > 2, there exists C Ξ,ς,ε > 0 such that III ≤ C Ξ,ς,ε G −1 = o(m −1 ). Also from Lemma E.3 , there exists C Ξ,ς,ε > 0 such that . Combining the bounds on III-V , the proof is complete.
Proof. Under Assumption 4.1 (i), Proposition 6.2 of Zhang and Wu (2021), combined with the arguments adopted in the proof of their Theorem 4.1 (most notably, their Equation (B.15)) and Bonferroni correction, obtains that there exist universal constants C 1 > 0 and C ν,α , C α > 0 that depend only on their subscripts, such that ν,α δ ν/2−1 z ν/2 + C 1 (4m + 1)GK χ exp − δz 2 C α mΨ 4 4,α thanks to Lemma E.1. Then, as in the proof of Lemma E.7, we can find C > 0 independent of i, i and show the first part of the claim by Lemma E.1. Similarly, under Assumption 4.1 (ii), Lemma E.1 and Theorem 6.3 of Zhang and Wu (2021) show that there exists a universal constant C 2 > 0 such that , which completes the proof.

E.4 Proof of Theorem 4.2
We provide a series of supporting results under the assumptions made in Theorem 4.2, leading to the proof of the claims. In what follows, we operate in M χ n,p . We defině Also, let δ χ,k = θ χ,k+1 − θ χ,k for 0 ≤ k ≤ K χ , and δ χ, Proposition E.9. (a) There exists a constant C > 0 such that Proof of (a). Under Assumption 4.2 (ii), applying Proposition 6.2 and Theorem 6.3 of Zhang and Wu (2021) with their (B.15), there exist universal constants C 1 , C 2 > 0 not dependent on 1 ≤ i, i ≤ p and constants C α , C ν,α > 0 that depend only on their subscripts, such that for any z > 0, Next, we bound the bias term We can show that I = O(m −1 ) by the arguments analogous to those adopted in bounding the RHS of (E.10). Also on M χ n,p , we have for fixed ≥ 0, by Lemma E.3, from which we obtain II = O(ρ n,p /G). In summary, we obtain E max k sup ω σ x, θ χ,k+1 ,ii (ω, δ χ,k ) − σ x,θ χ,k+1 ,ii (ω, δ χ,k ) 2 ≤ C ψ n ∨ 1 m ∨ ρ n,p G 2 for some constant C > 0 that does not depend on 1 ≤ i, i ≤ p, from which the conclusion follows.
Proposition E.14. Under the assumptions made in Theorem 4.2, we have on M χ n,p , Proof. By definition, we have where, from Lemma E.16, we have I = O p (θ n,p ) withθ n,p defined therein. Also, on M χ n,p , from Theorem 4.2. The conclusion follows by noting thatθ n,p ∨θ n,p = O(ϑ n,p ).
A consequence of Proposition E.14 is that P(E (2) n,p ) → 1, where with M as in (14).
Proof. We first note that solving (10) is equivalent to solving the problem column-wise, i.e.