Estimation and Inference of Time-Varying Auto-Covariance under Complex Trend: A Difference-based Approach

We propose a difference-based nonparametric methodology for the estimation and inference of the time-varying auto-covariance functions of a locally stationary time series when it is contaminated by a complex trend with both abrupt and smooth changes. Simultaneous confidence bands (SCB) with asymptotically correct coverage probabilities are constructed for the auto-covariance functions under complex trend. A simulation-assisted bootstrapping method is proposed for the practical construction of the SCB. Detailed simulation and a real data example round out our presentation.


Introduction
Our discussion begins with a heteroscedastic nonparametric regression model where Y i are the observations, µ is an unknown mean function, t i = i/n are the design points, i = 1, . . . , n, ε i are the errors with mean zero and variance 1, and V is the variance function. Historically, it has been assumed that the errors ε i are independent. Variance estimation in regression models with the unknown mean has traditionally been a rather important problem. Accurate variance estimation is required for the purpose of, for example, construction of confidence bands for the mean function, testing the goodness of fit of a model, and also in order to choose the amount of smoothing needed to estimate the mean function; see e.g. [29], [13], [15], and [18]. An extensive survey of the difference sequence approach to estimate the variance in the nonparametric regression setting when the variance function is only a constant can be found in [7]. The situation when the variance is not constant is more complicated. One of the first attempts to estimate the variance function in a regression model was made in [20] who proposed the basic idea of kernel smoothing of squared differences of observations. This idea has been further developed in [22]. [2] introduced a class of difference-based local polynomial regression-based estimators of the variance function V (x) and obtained optimal convergence rates for this class of estimators that are uniform over broad functional classes. [31] obtained the minimax rate of convergence for estimators of the variance function in the model (1) and characterized the effect of not knowing the mean function on the estimation of variance function in detail. Similar approach was used to construct a class of difference-based estimators in [3] when the covariate X ∈ R d for d > 1.
All of the above mentioned papers only considered the case where the data are independent. However, difference-based methods have also been used to estimate variance and/or autocovariance in nonparametric regression where the errors are generated by a stationary process. The pioneering approach here was probably that of [21] who proposed estimators based on the first-order differences to estimate (invertible) linear transformations of the variance-covariance matrix of stationary m-dependent errors. Here, by m-dependent errors we mean the errors generated by a stationary process whose autocovariance is equal to zero for any lag greater than some m > 0. [19] suggested second order differences to estimate the zero frequency of the spectral density of stationary processes with short-range dependence. In the case of autoregressive errors, [16] proposed √ n-root consistent and, under the assumption of normality of errors, efficient estimators of the autocovariance that are also based on differences of observations. Under certain mixing conditions, [26] proposed estimating the autocovariance function by applying difference-based estimators of the first order to the residuals of a kernel-based fit of the signal. [34] provided an optimal difference-based estimator of the variance for smooth nonparametric regression when the errors are correlated. Finally, the closest to us in spirit is, probably, [30] that proposed a class of difference-based estimators for the auto-covariance in nonparametric regression when the signal is discontinuous and the errors form a stationary m-dependent sequence. To the best of our knowledge, the problem of auto-covariance estimation in a nonparametric regression where the errors form a non-stationary sequence while the signal is discontinuous has not been considered before.
The purpose of this article is to estimate and make inference of the time-varying covariance structure of a locally stationary time series when it is contaminated by a complex trend function with both smooth and abrupt changes. Here local stationarity refers to the slowly or smoothly evolving data generating mechanism of a temporal system ( [6], [24], [37]). In time series analysis, the estimation and modelling of the auto-covariance structure is of fundamental importance in, for example, the optimal forecasting of the series ( [1]), the efficient estimation of time series regression models ( [17]) and the inference of time series regression parameters ( [1]). When the trend function is discontinuous, removing the trend from the time series and then estimating the auto-covariances from the residuals is not a good idea since it is very difficult to estimate the trend function near the points of discontinuity accurately. In this case, the aforementioned difference-based methods offer a good alternative. In this paper, we adopt a difference-based local linear regression method for the aforementioned time-varying auto-covariance estimation problem. The method can be viewed as a nonparametric and non-stationary extension to [30]. It is shown that the uniform convergence rate of auto-covariance function estimation for the difference-based method under complex trend is the same as that of auto-covariance function estimation of a zero-mean time series when the number of points of discontinuity as well as the jump sizes diverge to infinity at a sufficiently slow rate. Therefore, asymptotically, the accuracy of auto-covariance function estimation will not be affected by the complex trend when the difference-based nonparametric method is used.
Making inference of the auto-covariance functions is an important task in practice as practitioners and researchers frequently test whether certain parametric or semi-parametric models are adequate to characterize the time series covariance structure. For instance, one may be interested in testing whether the auto-covariance functions are constant over time so that a weakly stationary time series model is sufficient to forecast the future observations. There is a rich statistical literature on the inference of auto-covariance structure of locally stationary time series, particularly on the testing of weak stationarity of such series. See for instance [25], [12], [8], [23], [9] and [11]. To our knowledge, only constant or smoothly time-varying trend were considered in the aforementioned literature of covariance inference. In this paper, simultaneous confidence bands (SCB) with asymptotically correct coverage probabilities are constructed for the time-varying auto-covariance functions when estimated by the difference-based local linear method. The SCB serves as an asymptotically correct tool for various hypothesis testing problems of the auto-covariance structure under discontinuous mean functions. A general way to perform such hypothesis tests is to estimate the auto-covariance functions under the parametric or semi-parametric null hypothesis and then check whether the fitted functions can be fully embedded into the SCB. As the auto-covariance functions can be estimated with faster convergence rates under the parametric or semi-parametric null hypothesis, the aforementioned way to perform the test achieves correct Type-I error rate asymptotically. The tests are of asymptotic power 1 for local alternatives whose uniform distances from the null are of the order greater than that of the width of the SCB, see Theorem 2 in [35] for instance. We also propose a simulation-assisted bootstrapping method for the practical construction of the SCB.
The paper is organized as follows. In Section 2, we introduce the model formulation and some assumptions on x i and ε k i . Section 3 presents the asymptotic theory for local estimate β k (·). Practical implementation including a suitable difference lag and tuning parameters selection procedure, estimation of covariance matrices as well as an assisted bootstrapping method are discussed in Section 4. In Section 5, we conduct some simulation experiments on the performance of our SCBs. A real data application is provided in Section 6. The proofs of the main results are deferred to the Appendix.

Model formulation
Consider model: where µ i,n := µ(t i ) is a mean function or signal with unknown change points, y i,n = y(t i ) and x i,n = x(t i ) := G( i n , F i ) is a zero-mean locally stationary process with t i = i/n, i = 1, ..., n. Eq. (2) covers a wide range of nonstationary linear and nonlinear processes, see [37] for more discussion. We shall omit the subscript n in the sequel if no confusion arises. Let ζ i , i ∈ Z, be independent identically distributed (i.i.d.) random variables, and define F i = (..., ζ i−1 , ζ i ). Then, the process {x i } can be written as where G(·, ·) is a measurable function such that G(t, F i ) is well defined for all t ∈ [0, 1]. In this paper, we focus on the case that there exists 0 = a 0 < a 1 < · · · < a d < a d+1 = 1 such that where µ j (t) is a Lipschitz continuous function over [a j , a j+1 ) and d is the total number of change points. Till the end of this paper, we will always assume d = d n = O(n α ) and the maximal jump size ∆ n = max 1≤j≤d |µ j (a j ) − µ j (a − j )| = O(n β ) with 0 ≤ α, β < 1. To estimate the second order structure of the process Eq. (2), we introduce the approach based on the difference sequence of a finite order applied to the observations y i . Assuming that the number of observations is n + k, this difference-based covariance estimation approach would define simple squared differences of the observations, i.e., ρ k i := ρ k (t i−k ) = (y i − y i−k ) 2 , i = k + 1, k + 2, ..., k + n. Notice that for any fixed t, {ξ i (t) := G(t, F i )} i∈Z is a stationary process. For convenience, let us denote s i = t i−k = i−k n . Then, γ k (s i ) is the kth order autocovariance function of the process {x i } at the fixed time s i ; in other words, γ k (s i ) := Cov(x i , x i−k ), k = 0, 1, .... If k = 0, then γ 0 (s i ) turns out to be the variance of x i . We first introduce some notation that will be used throughout this paper. For any as the function space on [0,1] of functions that have continuous first p derivatives with integer p > 0. Now, we need the following definition and assumptions: Definition 1 (Physical dependence measure). Let (ζ j ) j∈Z be an i.i.d. copy of (ζ j ) j∈Z . Then, for any j ≥ 0, we denote F j = (F −1 , ζ 0 , ζ 1 , ..., ζ j ). The physical dependence measure for a stochastic system L(t, F j ) is defined as If j < 0, let δ q (L, j) = 0. Thus, δ q (L, j) measures the dependence of the output L(t, F j ) on the single input ζ 0 ; see [32] for more details.
Assumption 2 (Stochastic Lipschitz continuity). There exists a constant C > 0, such that G(t, Assumption 1 shows that the dependence measure of time series {x i } decays at a polynomial rate, thus indicating short-range dependence. Assumption 2 means that G changes smoothly over time and ensures local stationarity. Here, we show some examples of the locally stationary linear and nonlinear time series that satisfy these assumptions. Example 1 (Nonstationary linear processes). Let ζ i be i.i.d. random variables with ζ i ∈ L q , q ≥ 1; let a j (·), j = 0, 1, ..., be C 1 ([0, 1]) functions such that is well defined for all t ∈ [0, 1]. Clearly by [37, Proposition 2], we know that Assumption 1 will be satisfied if sup the stochastic Lipschitz continuity condition in Assumption 2 also holds true.
Example 2 (Nonstationary nonlinear processes). Let ζ i be i.i.d. random variables and consider the nonlinear time series framework where R is a measurable function and t ∈ [0, 1]. This form has been introduced by [37] and [36]. Suppose that for some x 0 , we have sup t∈[0,1] R(t, x 0 , ζ i ) q < ∞ for q > 0. Denote It is known from [37, Theorem 6] that if χ < 1, then Eq. (5) admits a unique locally stationary solution with ξ i (t) = G(t, F i ) and the physical dependence measure satisfies that δ q (G, j) ≤ Cχ j , which shows geometric moment contraction. Hence, the temporal dependence with exponentially decay indicates Assumption 1 holds with q = 8. Further by [37,Proposition 4], we conclude that Assumption 2 holds for q = 4 if Due to the local stationarity of the process {x i }, we have the following lemma which shows that, under mild assumptions, the auto-covariance of {x i } also exhibits polynomial decay.
With the above result, we can choose h large enough such that γ k (t) ≈ 0 for k ≥ h. Next we focus on the difference series ρ k i for k = 1, ..., h and we always assume h = h n = O(n 1/4 log n). By Eq. (2), we know that Recall s i = (i − k)/n for i = k + 1, ..., .k + n and notice that α k i := α k (s i ) = (x i − x i−k ) 2 is the squared difference of two locally stationary processes. Therefore, it is also a locally stationary process. As a result, we can define where β k (·) is the unknown trend function and ε k i := ε k (s i ) is a zero-mean process. Then ε k i can be written as where H k is a measurable function similar to G. With Eq. (7), if the trend function is smooth, one can easily obtain the estimator of β k (·). Now, we introduce the following conditions.
Assumption 3 guarantees that the trend function β k (·) changes smoothly for each k = 1, ..., h and is three-times continuously differentiable over [0, 1]. Assumption 4 prevents the asymptotic multicollinearity of regressors. Assumption 5 allows popular kernel functions such as Epanechnikov kernel. Now substituting Eq. (7) to Eq. (6), we have Since the length of the series {ρ k i } k+n i=k+1 is n, we reset the subscript with respect to i as {ρ k i } n i=1 and therefore the time point turns out to be t i = i/n for i = 1, ..., n. Similar notations are used for series {ε k i }, {λ k i } and {θ k i }. By Assumption 3 and the Taylor's expansion on β k (·), it is natural to estimate β k (t) using the local linear estimator as follows: where t i = i/n with i = 1, ..., n and K b (·) = K(·/b) is a kernel function, b = b n is the bandwidth satisfying b → 0 and nb → ∞. Since Eq. (11) is essentially a weighted least squares estimate, we can write the solution of Eq. (11) as where is the weight given to each observation.
Next, we will establish the following two lemmas that are useful in establishing asymptotic properties of proposed estimators. Their proofs are given in the Appendix.

Asymptotic theory
By Assumption 3 and for l = 0, 1, ..., define Then Eq. (11) can be expressed as Let Theorem 1. Suppose that Assumptions 1-5 hold and further assume that Let us comment on the conditions listed in Theorem 1. Condition (1) shows the smoothness of σ k (t). Condition (2) indicates that the change-point number and size can both go to infinity but at a slow rate. The assumption nb 5 log(n) → 0 in Condition (3) is an undersmoothing requirement that reduces the bias of the estimators to the second order. Notice when there exists at least a change point on µ(·). However, the estimate of λ k i can be viewed as a negligible term (see Eq. (21) in the proof of Theorem 2). With the previous discussion in mind, we can define where b h and b k are the bandwidths for estimators γ 0 (t) and γ k (t), respectively. Making it easy to distinguish, here we use the different notations for the bandwidths which will be selected by some criterion (see Section 4.4). Notice that we require the same bandwidth (b k ) to compute the estimator of γ k (t). With the above results, the SCB for γ 0 (·) is straightforward.
Corollary 2. Suppose that the smallest eigenvalue of σ C,k (t) is bounded away from 0 on [0, 1] for k = 1, ..., h − 1. Moreover, we assume that all of the conditions of Theorem 1 are valid. Then, we have (i) Remark 1. It is noteworthy to mention that for estimating β C,k (t), we use the same bandwidth b k ; therefore, the entire estimator β C, depends on only a single tuning parameter (bandwidth b k ). This enables us to achieve the conclusion of Corollary 2(i) based on the result of Theorem 1. As a result, Corollary 2(ii) also holds true due to this fact.
After constructing SCBs for the second-order structure γ k (·), the following theorem states that γ k (t) are consistent estimators for γ k (t) uniformly in t for all k = 0, ..., h − 1.

Theorem 2. Under Assumptions 1-5 and suppose conditions
This theorem implies the uniform consistency of β k,b (·). Additionally, due to the relationship between β k and γ k , we can also easily obtain the following consistency result for γ k (·).
Corollary 3. With the conditions in Theorem 2, we have 4 Practical implementation 4

.1 Selection of the difference lag
Note that for any fixed time t, β k,b k (t) = 2 γ 0 (t) − 2 γ k (t) and recall that γ k (t) ≈ 0 when k ≥ h, where h is a large value that has been chosen in advance. Hence, we know that if k ≥ h, β k,b k (t) ≈ 2 γ 0 (t) is practically invariant with respect to t as k increases. This fact suggests the following bandwidth selection procedure.
First, for any fixed t, we choose a large enough value h 0 (t) and select k = h 0 (t). Next, we calculate β k,b k (t). Then, by successively decreasing the value of k and considering k = h 0 (t) − 1, h 0 (t) − 2, ..., we calculate the corresponding quantities β k,b k (t) until β k,b k (t) shows an abrupt change. At this point, the optimal difference lag for time t can be selected as the current k plus 1. Intuitively, we can interpret this through the scatterplot of (k, β k,b k (t)). When the slope of the function β k,b k (t) shows an obvious change, then we can choose h * (t) = k + 1. Following the above procedure for each time point t i = i/n, i = 1, ..., n, we finally choose the optimal lag as h = n i=1 h * (t i )/n.

Covariance matrix estimation
To apply Corollaries 1 and 2 (ii), we need to estimate the long-run variance Σ 2 k (·) in Eq. (16) first. This problem is complicated but has been extensively studied by many researchers. Here we adopt the technique considered by [38].
Let Q k i = m j=−mε i+j , i = 1, ..., n, whereε i = (ε h i , ε k i ) for i = 1, ..., n. Notice that E(ε i ) = 0 and denote N k i := Q k i Q k i /(2m + 1). In the locally stationary case, we can make use of the fact that a block of {ε i } is approximately stationary when its length is small compared with n. Hence, E(N k i ) ≈ Σ 2 k (t i ) as m → ∞ and m/n → 0. Let τ = τ n be the bandwidth and define the covariance matrix estimator as , with τ being the bandwidth. Therefore, the estimate Σ 2 k (t) is guaranteed to be positive semidefinite. The following theorems provide consistency of our covariance matrix estimate.
, m = m n → ∞, m = O(n 1/3 ), τ → 0 and nτ → ∞. Then, (i) for each k, k = 1, ..., h − 1 and any fixed t ∈ (0, 1), In practice, the errorsε j cannot be observed, thus we use Σ 2 i is defined as N k i withε j therein replaced by its estimatorε j . Theorem 4. Assume that conditions of Theorem 2 and conditions of Theorem 3 hold. Denote ν n = m log(n)χ n , where χ n is defined in Theorem 2 and further assume ν n → 0. Then Note that σ 2 h (t) is the first diagonal element of Σ 2 k (t) and σ 2 C,k (t) = C Σ 2 k (t)C. Thus, the covariance estimates in Corollaries 1 and 2 (ii) can be easily calculated via plugging in the long-run covariance matrix estimate Σ 2 k (t).

Simulation assisted bootstrapping method
Now we aim to apply Corollary 1 and Corollary 2 (ii) to construct the SCBs. Let γ 0 (t) and γ k (t) be uniformly consistent estimators of γ 0 (t) and γ k (t) for k = 1, ..., h−1, respectively. Then the corresponding (1 − α)th SCB with α ∈ (0, 1) for γ 0 (t) and γ k (t) are Due to the slow rate of convergence to Gumbel distribution, in practice, the UCB from Corollaries 1 and 2 (ii) may not have good finite-sample performances. To circumvent this problem, we shall adopt a simulation assisted bootstrapping approach. Proposition 1. Suppose conditions in Theorem 1 hold and also assume that σ k (t) is Lipschitz continuous for k = 1, ..., h. Then, on a richer probability space, there are i.i.d. standard normal distributed random variables u i such that , which can be obtained by generating a large number of i.i.d. copies via bootstrapping. Therefore, the above proposition provides us with an alternative way to construct the SCB of the autocovariance function without using the asymptotic Gumbel distribution.
For ease of application, we combine procedures mentioned above into a convenient sequence of steps below.
• Choose the difference lag order h by using method that is proposed in Section 4.1.
• Generate i.i.d. random variables u 1 , u 2 , ... ∼ N (0, 1) and calculate sup t∈ • Repeat the last step for a large number of times (e.g. 10 4 ) and obtain the estimated • Calculate Σ 2 k (t) by using the method in Section 4.2. Then, obtaining σ h (t) together with σ C,k (t) is straightforward.

Selection of tuning parameters
In this subsection, we briefly discuss the practical choices of tuning parameters b, m and τ . Here, we consider the generalized cross-validation (GCV) method by [5] to choose the bandwidth b. Specifically, we consider two cases of bandwidth selection for γ 0 (t) and γ k (t), k = 1, ..., h − 1, respectively. For estimating γ 0 (t), let P h = (ρ h 1 , ..., ρ h n ) and P h (b) = (ρ h 1 (b), ...,ρ h n (b)) be the corresponding fitted values. One can write P h (b) = H(b)P h , where H(b) is an n by n square hat matrix that depends on b. Then, we choose the optimal bandwidth (say b h ) that minimizes On the other hand, when estimating γ k (t) for k = 1, ..., h − 1, we treat γ k (t) as a whole term and choose a joint bandwidth for it. Similarly, denote P k = (ρ h 1 −ρ k 1 , ..., ρ h n −ρ k n ) and ) be the corresponding fitted values. As before, one can write P k (b) = H(b)P k . With this in mind, we select as optimal the bandwidth (say b k ) that minimizes the following quantity: For the choice of m and τ , we now employ the extended minimum volatility method (including two parameters) which was proposed in [10,Chapter 9]. This method is based on the fact that if a pair of block size and bandwidth is in an appropriate range, then confidence regions for the local mean constructed by Σ 2 k (t) should be stable. Therefore, we first consider a grid of possible block sizes and bandwidths and then choose the optimal pair that minimizes the volatility of the boundary points of the confidence regions in the neighborhood of this pair. To be more specific, let the grid of possible block sizes and bandwidths be {m 1 , ..., m M 1 } and {τ 1 , ..., τ M 2 }, respectively. Then denote the estimated long-run covariance matrices as { Σ 2 k (m i , τ j , t)} for i = 1, ..., M 1 , j = 1, ..., M 2 . For each pair (m i , τ j ), we need to calculate where ISE denotes the integrated standard error

Simulations
To illustrate performance of the proposed estimator of autocovariance, we consider several models. For each model, we obtain the uniform confidence interval coverage of the true variance function and the autocovariance function at lag 1 for three different sample sizes: n = 400, n = 600 and n = 800. In each case, we use 500 replications. To select bandwidths b h and b k we use the grid from 0.15 to 0.45 with the step size 0.01. We also provide a graphical illustration of a confidence interval enclosing the true variance and autocovariance lag 1 functions for each of the models considered. The first model considered has the errors that are generated by a locally stationary linear process (4) with a j (t) = t 2 j , j = 1, 2, . . . while the sequence (ζ i ) consists of iid normal random variables with mean zero and variance 1. In this case the coefficients start with j = 1 since otherwise a 0 (t) is undefined at t = 0. The Assumption 1 is satisfied since sup t∈[0,1] a 2 j (t) ≤ 1 4 j which is, of course, O(j −2 ). Assumption 2 is also satisfied because ∞ j=1 sup t∈[0,1] a 2 j (t) = ∞ j=1 j 2 4 j < ∞. Next, the mean function µ(t) is taken to be a piecewise constant function with six change-points located at fractions 1 6 ± 1 36 , 3 36 ± 2 36 , and 5 6 ± 3 36 of the sample size n. In the first segment, µ 1 ≡ µ(t 1 ) = 0, in the second it is equal to 1, and in the remaining segments µ(t) alternates between 0 and 1, starting with 0 in the third segment. This mean function is very similar to the one that has been considered earlier in several other publications; see e.g. [4] and [30].
The second model we consider has exactly the same error structure as Model 1 but the mean function is a slightly different one. In particular, we make the value of the function in the second segment 2 instead of 1 while the remaining segments of µ(t) alternate between 0 and 1, starting with 0 in the third segment. Since the error process remains the same as before in Model 1, Assumptions 1 and 2 are satisfied.
The third model we consider is where the errors are generated by a locally stationary MA(2) process a j (t)ζ i−j with coefficients being equal to a j (t) = (t+0.05) j 2 j , j = 0, 1, 2. The sequence (ζ i ) consists of iid N (0, 0.3) random variables. The locally stationary MA process considered is a special case of the general locally stationary linear process. Since the process consists of the finite number of terms, the stochastic Lipschitz continuity condition in the Assumption 2 is satisfied automatically. Because sup t∈[0,1] a 2 j (t) ≤ (0.525) 2j , the Assumption 1 will also be satisfied. Finally, the mean function stays the same as in the Model 1.
In Tables 1 to 3, we illustrate coverage probabilities of uniform confidence intervals of the variance function and lag 1 covariance function for all three of the models considered. We also consider three possible sample sizes, n = 400, n = 600 and n = 800. Note that even a relatively small sample size of 400 gives excellent coverage probabilities. It is also worthwhile noting that the coverage probabilities are generally higher for lag 1 autocovariance function than for the variance function.   Table 3: Empirical coverage probabilities for all the models when the sample size is 800 To illustrate the behavior of uniform confidence intervals for each of the three models considered, we also include sample plots of fitted variance/autocovariance curves with corresponding confidence intervals. For each model, two plots are given: one with the true variance function, its estimate, and a uniform confidence interval for the estimated variance curve, while the other one contains the true lag 1 autocovariance function, its estimate, and the corresponding uniform confidence interval for the estimated autocovariance curve. In each of the plots, a solid line is used for the true variance/autocovariance curve, a dashed line for the corresponding estimated curve, and red dotted lines for uniform confidence intervals.
To illustrate the reasonableness of our method, we also provide a quick comparison of our approach to a very straightforward "naive" method. Such a method would start with a rough estimate of the mean function µ(t) using a local smoother, for example, a local linear regression. The resulting rough estimate of the mean function can then be subtracted from observations to form a series of residuals e i , i = 1, . . . , n. Using this series, a naive approach would estimate the variance function γ 0 (t) by applying a smoother, e.g. yet again the local linear regression, to squared residuals e 2 i , i = 1, . . . , n. In much the same way, applying the local linear regression to a series e i e i−1 , i = 2, . . . , n will result in a "naive" estimate of the lag 1 autocovariance function γ 1 (t). In both situations, we used a simple generalized cross-validation to obtain the optimal smoothing bandwidth.
It is probably sufficient to say that such a naive approach fails completely in an attempt to estimate the second order structure when the mean is discontinuous and has numerous change points. More specifically, we tried to obtain the coverage of the true variance function γ 0 (t) by a uniform confidence interval that is based on the "naive" estimate described above. To do so, we used our Model 1 with the sample size n = 400 and 500 replications. We found that the coverage is zero, that is, the true variance function γ 0 (t) is never completely inside the uniform confidence interval. This can be explained properly by noticing that our mean estimate used to obtain residuals is extremely crude. More specifically, in order for the local linear mean function estimator to be consistent at a given point the mean function µ(t) has to have two continuous derivatives at that point; see e.g. [14] p. 62 for a detailed discussion. This lack of consistency results in a severe bias of the variance function estimator. Thus, such a direct approach seems to be completely inappropriate for determination of the second order structure.

Real data application
In this section, we illustrate our approach using a real dataset. There is a rather clear evidence that the global temperatures are nonstationary (see e.g. [28]) and so we use the dataset that consists of monthly temperature anomalies observed during the period from January 1856 to September 2019. A shorter subset of the same series has been used earlier in [27]. The data used are publicly available from the Climate Research Unit of the University of East Anglia, UK at https://www.cru.uea.ac.uk/. The anomalies are defined here as the difference of temperatures from a reference value. The anomaly data are available for both Northern and Southern hemisphere separately. Figures (4a) and (4b) display the temperature anomaly data for both hemispheres.
Our purpose is to estimate the variance and lag 1 autocovariance function of these data as a function of time. For the Northern hemisphere data the approach suggested in our manuscript produces an almost monotonically decaying variance curve that suggests that some nonstationarity is, indeed, present in the data. This monotonic decay is probably due to the increasing number of weather stations recording the data over time. The variance of the Southern hemisphere data is also mostly decreasing although the decay is not as clearly monotonic as for the Northern hemisphere data. Note that, for both sets of data, the lag 1 autocovariance is very small in magnitude; however, the horizontal zero  line added to both autocovariance plots is clearly not fully inside the uniform confidence interval, indicating that the temperature series are not white noises. Now, we are interested in testing whether for these data there exist change points in the mean function, namely, H 0 : µ 1 = µ 2 = · · · = µ n = µ ←→ H 1 : µ i = µ j for some 1 ≤ i < j ≤ n.
To this end, we will use the robust bootstrap test for nonstationary time series proposed by [36][Section 4]. For Northern Hemisphere data, the robust bootstrap test yields a < 0.1% p-values with 10000 bootstrap samples, which provides a very strong evidence against the null hypothesis of no structural change in mean. On the other hand, we applied the robust bootstrap to the Southern Hemisphere data. The corresponding p-value of the test with 10000 bootstrap samples is also < 0.1%, which also shows a strong evidence against As a result, the test further illustrates the usefulness of our method for constructing SCB with finite change points. Over some time periods, the data with wild fluctuations indicates a change in mean and suggests the non-stationarity, as pointed by [28].

Appendix
The following theorem provides the Gaussian approximation result for nonstationary multiple time series, which can be found in Theorem 2 of [38].
Theorem 5 (Theorem 2 in [38]). Suppose that Assumptions 1 and 2 hold. A partial sum process can be defined asS k i = i j=1 ε k j for k = 1, ..., h. Then, on a richer probability space, there exist i.i.d. standard normal random variables u 1 , u 2 , ..., and a process S k i such where σ k (·) is defined as Eq. (9).
Theorem 5 implies the Gaussian approximation for a partial sum of a locally stationary process. Note that, due to the result stated in Lemma 2, the physical dependence measure δ 4 (H k , i) has different types of the polynomial decay under two circumstances, which is more complicated than that of Corollary 2 in [33]. But letting the order of the m-dependence sequence larger than k and making a careful check of the proof of Corollary 2 in [33], we can obtain the same argument. Owing to the non-stationarity, the approximated Gaussian process { i j=1 σ k (t j )u j } n i=1 has independent but possibly non-identically distributed increments.
Proof of Lemma 1. For j ∈ Z, define the projection operator P j (·) = E(·|F j ) − E(·|F j−1 ), then we can write x i = i j=−∞ P j (x i ). Denote t i = i/n, we have The first inequality follows by the orthogonality of P j (·) and the second inequality is due to Fubini's theorem and Cauchy-Schwartz inequality. The last inequality follows from the argument in [32,Theorem 1]. Therefore, with Assumption 1, there exists a constant C > 0 such that |γ k (t i )| ≤ Ck −2 . This concludes our proof.
Proof of Lemma 2. Now, we consider the locally stationary process (x l − x l−k ) 2 and let x l be the coupled process of x l with ζ 0 replaced by an i.i.d. copy ζ 0 . Then for each k = 1, ..., h, Proof of Lemma 3. By Assumption 3, we know that β k (·) is also Lipschitz continuous, thus with Eq. (7), ≤C|t − s|.
The first inequality follows from the triangle inequality and Minkowski's inequality.
To prove Theorem 1, we need to introduce the following lemmas.
Proof. Similar to Lemma 2 in [38]. Next, we treat T k n,0 (t) as D k (t) in Lemma 5, then Theorem 1 follows.
Proof of Theorem 2. Recall that our model contains d = O(n α ) change points with the maximal size ∆ n = O(n β ) on µ(·) and let Ω and Λ k be the n-dimensional vectors with the entrywises ω b n (t, i) and (µ i − µ i−k ) 2 for i = 1, ..., n, respectively. Due to the fact that µ(·) is Lipschitz continuous, one can see that Λ k consists of kd components being O(n 2β ) and other components being O(1/n 2 ). For each k = 1, ..., h and any fixed t ∈ [0, 1], by Eqs. (6), (7) and (12) we have ω b n (t, i)θ k i := I + II + III + IV.
It is obvious that I = O P (b 2 ). Then we will apply a chaining argument for calculating II. As for III, notice that In the end, by the similar chaining argument as those in the proof of I, we have Note that the assumption α + 2β < 2/5 entails β < 1/5, therefore by elementary calculation, the above four kinds of bounds all converge to 0 as n → ∞.
Proof of Theorem 3. Similar to the proof of Theorem 4 in [38].
Proof of Theorem 4. Let I be a closed interval in (0, 1) such that I ⊂ I and the two intervals do not share common end points. Recall β k (t) = (β h (t), β k (t)) and denote β k,b k (t) = ( β h,b k (t), β k,b k (t)) . According to Theorem 2, it follows that Note that Q k i /(2m + 1) is the Nadaraya-Waston smoother of the series {ε i } at i with the rectangle kernel and bandwidth m/n. Therefore, for each k = 1, ..., h, we have sup i/n∈I Let Q k i = m j=−mε i+j and N k i = Q k i Q k i /(2m + 1). Then Substituting equations (22) and (23) into the above equation, we have sup i/n∈I |N k i − N k i | = O P (ν n ) with the assumption ν n → 0. By the definitions of Σ k (t) and Σ k (t), we obtain sup Together with the results of Theorem 3, Theorem 4 holds.
Proof of Proposition 1. This proposition follows by Theorem 1 and Eq. (19) from the proof of Lemma 5.