\renewtheoremstyle

plain ProofProof \nolinenumbers \jyear \jvol \jnum

On the partial autocorrelation function for locally stationary time series: characterization, estimation and inference

Xiucai Ding xcading@ucdavis.edu Department of Statistics, University of California, Davis, Davis 95616 USA Zhou Zhou zhou.zhou@utoronto.ca Department of Statistical Sciences, University of Toronto, Toronto M5G 1X6, Canada

Abstract

For stationary time series, it is common to use the plots of partial autocorrelation function (PACF) or PACF-based tests to explore the temporal dependence structure of such processes. To our best knowledge, such analogs for non-stationary time series have not been fully established yet. In this paper, we fill this gap for locally stationary time series with short-range dependence. First, we characterize the PACF locally in the time domain and show that the $j$ th PACF, denoted as $\rho_{j}(t),$ decays with $j$ whose rate is adaptive to the temporal dependence of the time series $\{x_{i,n}\}$ . Second, at time $i,$ we justify that the PACF $\rho_{j}(i/n)$ can be efficiently approximated by the best linear prediction coefficients via the Yule-Walker’s equations. This allows us to study the PACF via ordinary least squares (OLS) locally. Third, we show that the PACF is smooth in time for locally stationary time series. We use the sieve method with OLS to estimate $\rho_{j}(\cdot)$ and construct some statistics to test the PACFs and infer the structures of the time series. These tests generalize and modify those used for stationary time series in Brockwell & Davis (1987). Finally, a multiplier bootstrap algorithm is proposed for practical implementation and an $\mathtt{R}$ package $\mathtt{Sie2nts}$ is provided to implement our algorithm. Numerical simulations and real data analysis also confirm usefulness of our results.

keywords:

Locally stationary time series; PACF; Sieve method; Multiplier bootstrapping.

\arabicsection Introduction

The partial autocorrelation function (PACF) is one of the most popular and powerful tools for stationary time series modelling and analysis Brockwell & Davis (1987). However, in the era of big data, as increasingly longer time series are being collected, it has become more appropriate to model many of those series as locally stationary processes whose data generating mechanisms evolve smoothly over time. In this setting, the effectiveness of the classical PACF deteriorates and it is of urgent demand to establish the theories of PACF for locally stationary time series.

Even though there exists a rich body of literature on locally stationary time series analysis, see Dahlhaus (2012); Dahlhaus et al. (2019) for a review, much less has been studied related to the PACF. In Dégerine & Lambert-Lacroix (2003), the authors generalized the characterization of PACFs via some useful decomposition as introduced in Ramsey (1974) to general non-stationary processes. They also briefly discussed how to estimate the PACFs based on a generalized Levinson-Durbin algorithm when the autocovariance function is given. More recently, in Killick et al. (2020), by generalizing the partial autocorrelations of stationary processes to locally stationary time series from a wavelet spectrum perspective, the authors provided two new estimators for the local PACFs. The consistency of the wavelet-based estimator and the asymptotic distribution of the windowed estimator under Gaussian assumption have also been studied. However, the decay speed of local PACFs as a function of the lags has not been established and a direct time-domain characterization of the PACFs of locally stationary processes has not been fully investigated. Moreover, the inference for PACFs of locally stationary time series, for example significance tests and PACF-based Portmanteau tests, are still missing in the literature.

Motivated by the above challenges, in this paper, we aim to systematically study the theories of PACFs for locally stationary time series. For characterization, in contrast to Dégerine & Lambert-Lacroix (2003); Killick et al. (2020), we define PACFs for general locally stationary time series in the time domain using stationary approximations at each time point (c.f. Definition \arabicsection.\arabicthm). There are several advantages in using this characterization. First, since the time series is approximately stationary locally, the lower-order PACFs can be well approximated by the best short-term linear prediction coefficients via the Yule-Walker equations; see (\arabicsection.\arabicequation) and (\arabicsection.\arabicequation). This connection not only allows us to study the PACF via ordinary least squares (OLS) estimation but also enables us to establish the decay properties of the PACFs which are adaptive to the temporal dependence decay of the time series; see (\arabicsection.\arabicequation). Second, the smoothness of the locally stationary covariance structure can be easily translated to that of the PACFs; see (\arabicsection.\arabicequation). Therefore, it suffices to estimate some smooth functions $\rho_{j}(\cdot)$ at different lags $j.$ More specifically, together with the OLS, the smooth PACFs can be estimated adaptively using the nonparametric method of sieves via flexible choices of basis functions such as the wavelets and the orthogonal polynomials Chen (2007); see Section \arabicsection.\arabicsubsection. Theoretically, under mild assumptions, the estimators are consistent uniformly in the time domain (c.f. Theorem \arabicsection.\arabicthm). Third, based on the time-domain characterization and the OLS form of the sieve estimators, one can further conduct various tests on the PACFs. For example, one can perform a white noise Portmanteau test or significance tests on some PACFs (e.g., checking the order for an AR process) uniformly over time and lags. Both tests have not been fully studied yet under the locally stationary time series framework. We establish the asymptotic normality and conduct a power analysis for the tests (c.f. Theorems \arabicsection.\arabicthm and \arabicsection.\arabicthm). We also propose a multiplier bootstrap procedure for practical implementation (c.f. Algorithm 1). Numerical simulations and real data analysis are provided to illustrate the usefulness of our results and an $\mathtt{R}$ package $\mathtt{Sie2nts}$ is provided for users. Since our method covers stationary time series as a special case, we promote the use of the $\mathtt{Sie2nts}$ package instead of the default $\mathtt{pacf}$ function in $\mathtt{R}$ which only handles stationary time series.

The paper is organized as follows. In Section \arabicsection, we provide the characterization of the PACFs for locally stationary time series and study their asymptotic properties. In Section \arabicsection, we introduce our estimator for PACFs based on the nonparametric sieve method and inference procedures based on multiplier bootstrap. In Section \arabicsection, we provide theoretical analysis for our estimation and inference procedures. Numerical simulations and real data analysis are offered in Section \arabicsection. More details are provided in our online supplement Ding & Zhou (2024). Especially, technical proofs are deferred to Section A, further discussions and remarks are provided in Section B, some tuning parameter selection algorithm is listed in Section C and additional simulation results are offered in Section D.

Conventions. For a random variable $x\in\mathbb{R}$ and some constant $q\geq 1,$ we denote by $\|x\|_{q}:=(\mathbb{E}|x|^{q})^{1/q}$ its $L^{q}$ norm. We simply write $\|x\|\equiv\|x\|_{2}$ when $q=2.$ For any deterministic vector $\bm{x}=(x_{1},x_{2},\cdots,x_{p})^{*}\in\mathbb{R}^{p},$ we use $|x|=\sqrt{\sum_{i=1}^{p}x_{i}^{2}}$ for its $\ell_{2}$ (or Euclidean) norm. For any matrix $A,$ we use $\|A\|$ to stand for its operator norm. For two sequences of deterministic positive values $\{a_{n}\}$ and $\{b_{n}\}$ , we write $a_{n}=\mathrm{O}(b_{n})$ if $a_{n}\leq Cb_{n}$ for some positive constant $C>0$ . Moreover, we write $a_{n}=\mathrm{o}(b_{n})$ if $a_{n}\leq c_{n}b_{n}$ for some positive sequence $c_{n}\downarrow 0.$ For a sequence of random variables $\{x_{n}\}$ and positive real values $\{a_{n}\},$ we use the notation $x_{n}=\mathrm{O}_{\ell^{q}}(a_{n})$ to state that $x_{n}/a_{n}$ is bounded in $L^{q}$ norm; that is $\|x_{n}/a_{n}\|_{q}\leq C$ for some finite constant $C>0.$ If $q=1,$ we simply write $x_{n}=\mathrm{O}_{\mathbb{P}}(a_{n}).$ We use $\mathbf{I}_{c}$ for a $c\times c$ identity matrix. We use $C^{d}([0,1])$ for the function space on $[0,1]$ of continuous functions that have continuous first $d$ derivatives.

\arabicsection Characterization of PACF for locally stationary time series

In this section, we provide the characterization of PACF of locally stationary time series and study its properties. Suppose that we observe $\{x_{i,n}\}_{i=1}^{n}$ . For simplicity and without loss of generality, we assume the time series is centered (i.e., mean is zero). Till the end of the paper, for notional simplicity, we always write $x_{i}\equiv x_{i,n}.$ From time to time, we will emphasize the dependence on $n$ for various quantities.

In this paper, we focus on a general class of locally stationary time series following definition introduced in Ding & Zhou (2023). It covers many commonly used locally stationary time series models in the literature. See for instance Ding & Zhou (2020); Dahlhaus et al. (2019); Dahlhaus (2012); Roueff & Sanchez-Perez (2018); Kley et al. (2019); Dette et al. (2011); Dette & Wu (2020); Vogt (2012); Zhou & Wu (2009). We refer the readers to Example B.\arabicthm in our supplement for more details.

Definition \arabicsection.\arabicthm (Locally stationary time series and its PACF).

A non-stationary time series $\{x_{i}\}$ is a locally stationary time series (in covariance) if there exists a function $\gamma(t,k):[0,1]\times\mathbb{N}\rightarrow\mathbb{R}$ such that

\operatorname{Cov}(x_{i},x_{j})=\gamma(t_{i},|i-j|)+\mathrm{O}\left(\frac{|i-j% |+1}{n}\right),\ t_{i}=\frac{i}{n}.

(\arabicsection.\arabicequation)

Moreover, we assume that $\gamma$ is Lipschitz continuous in $t$ and for any fixed $t\in[0,1],$ $\gamma(t,\cdot)$ is the autocovariance function (ACF) of some stationary process whose $j$ th order PACF is denoted as $\rho_{j}(t)$ . We shall also call $\rho_{j}(t)$ the $j$ th order PACF at rescaled time $t$ of $\{x_{i}\}.$

Observe that Definition \arabicsection.\arabicthm essentially means that the covariance structure of $\{x_{i}\}$ can be well approximated locally by that of a stationary process. Consequently, at each rescaled time $t$ , the PACF of $\{x_{i}\}$ is defined by that of the approximating stationary process. Now we study the PACF introduced above. Since only one realization of the series $\{x_{i}\}$ is available, we will need to assume short-range dependence that for some $\tau>1,$

\max_{k,n}\left|\operatorname{Cov}(x_{k,n},x_{k+r,n})\right|=\mathrm{O}\left(|% r|^{-\tau}\right)\ \text{and}\ \sup_{t}|\gamma(t,|r|)|=\mathrm{O}(|r|^{-\tau}),

(\arabicsection.\arabicequation)

and add additional regularity conditions. These will be summarized in Assumption \arabicsection.\arabicthm after some necessary notations are introduced. For the autocovariance function $\gamma(t,\cdot)$ in (\arabicsection.\arabicequation), given a lag $j,$ we define a vector of functions $\bm{\phi}_{j}(t)=(\bm{\phi}_{j,1}(t),\cdots,\bm{\phi}_{j,j}(t))^{*}:[0,1]% \rightarrow\mathbb{R}^{j}$ via the local Yule-Walker equation

\bm{\phi}_{j}(t)=\Gamma_{j}(t)^{-1}\bm{\nu}_{j}(t),

(\arabicsection.\arabicequation)

where $\Gamma_{j}(t)\in\mathbb{R}^{j\times j}$ is a symmetric matrix whose $(k,l)$ entry is defined as $\Gamma_{j}^{(k,l)}(t):=\gamma(t,|k-l|),$ and $\bm{\nu}_{j}(t)\in\mathbb{R}^{j}$ is vector whose $k$ th entry is defined as $\bm{\nu}_{j}^{(k)}(t):=\gamma(t,k).$ Consequently, we can write Brockwell & Davis (1987)

\rho_{j}(t)=\bm{\phi}_{j,j}(t).

(\arabicsection.\arabicequation)

Note that for stationary time series the PACF is closely related to the best linear forecast coefficients of the process (Brockwell & Davis (1987)). The representation (\arabicsection.\arabicequation) together with the smoothness of the time series covariance structure imply that we can study the PACF of locally stationary time series via local best linear forecasts and therefore a simple local regression analysis. Specifically, for all $1\leq j\leq i-1,$ denote the $j$ th order best linear forecast of $x_{i}$ as $\widehat{x}_{i,j}=\sum_{k=1}^{j}\phi_{ik,j}x_{i-k},$ where $\phi_{ik,j}\equiv\phi_{ik,j,n},\ 1\leq k\leq j,$ are the best linear forecast coefficients. Define the residual as $\epsilon_{i,j}\equiv\epsilon_{i,j,n}:=x_{i}-\widehat{x}_{i,j}.$ We now write

x_{i}=\sum_{k=1}^{j}\phi_{ik,j}x_{i-k}+\epsilon_{i,j}.

(\arabicsection.\arabicequation)

Using (\arabicsection.\arabicequation), we introduce the following notations

\rho_{i,j}=\phi_{ij,j},\ 1\leq j\leq i-1.

(\arabicsection.\arabicequation)

The following theorem studies the uniform decay property of the PACF and builds the connection between the PACF defined in Definition \arabicsection.\arabicthm and $\rho_{i,j}$ in (\arabicsection.\arabicequation) for the locally stationary time series.

Theorem \arabicsection.\arabicthm.

For the locally stationary time series $\{x_{i}\}$ satisfying Definition \arabicsection.\arabicthm, suppose Assumption \arabicsection.\arabicthm holds. Then we have the followings holds.

\arabicenumi.

For $\tau$ in (\arabicsection.\arabicequation), we have that for $\rho_{j}(t)$ in Definition \arabicsection.\arabicthm and $\rho_{i,j}$ in (\arabicsection.\arabicequation)

\sup_{i>j}|\rho_{i,j}|=\mathrm{O}\left(\left(\frac{\log j+1}{j}\right)^{\tau-1% }\right)\ \text{and}\ \sup_{t}|\rho_{j}(t)|=\mathrm{O}\left(\left(\frac{\log j% +1}{j}\right)^{\tau-1}\right).

(\arabicsection.\arabicequation)

\arabicenumi.

$\rho_{j}(t)\in C^{d}([0,1]),$ for some integer $d>0$ defined in (\arabicsection.\arabicequation) and

\sup_{i>j}\left|\rho_{j}\left(\frac{i}{n}\right)-\rho_{i,j}\right|=\mathrm{O}% \left(\frac{j^{2}}{n}\right).

(\arabicsection.\arabicequation)

Remark \arabicsection.\arabicthm.

Several remarks on Theorem \arabicsection.\arabicthm are in order. First, (\arabicsection.\arabicequation) implies that, uniformly over time, the PACF decays polynomially fast to 0 as a function of the lag with the speed adaptive to that of the autocovariance. If $\tau$ is sufficiently large, this implies that practically we only need to consider the first few lags of the PACF. In practice, people are usually concerned with the cutoff that $\sup_{t}|\rho_{j}(t)|\gg n^{-1/2}.$ In this sense, we conclude from (\arabicsection.\arabicequation) that we only need to focus on the lags for $j=\mathrm{O}\left(n^{\frac{1}{2(\tau-1)}}\right).$ Second, (\arabicsection.\arabicequation) demonstrates that under Definition \arabicsection.\arabicthm and Assumption \arabicsection.\arabicthm, the PACFs can be well approximated by $\rho_{i,j}$ if $j$ is not too large. Observe that $\rho_{i,j}$ is defined by local best linear forecasts and therefore is closely related to OLS. Therefore, (\arabicsection.\arabicequation) lays the theoretical foundation for the estimation and inference procedures in Section \arabicsection.

\arabicsection Estimation and inference procedures

In this section, we provide the procedures of the estimation and inference of the PACFs of locally stationary time series. The theoretical justifications will be provided in Section \arabicsection. As mentioned in Remark \arabicsection.\arabicthm, in what follows, we shall only consider lag $j\leq j_{*}$ with

j_{*}\asymp n^{\frac{1}{2(\tau-1)}}.

(\arabicsection.\arabicequation)

\arabicsection.\arabicsubsection Sieve nonparametric estimation

In this section, we estimate the PACFs $\rho_{j}(\cdot).$ As proved in Theorem \arabicsection.\arabicthm, since $\rho_{j}(t)\in C^{d}([0,1]),$ it is natural for us to approximate it via a finite diverging term basis expansion using the method of sieves as in Chen (2007); Ding & Zhou (2020, 2023). Recall $\rho_{j}(t)$ is defined via $\bm{\phi}_{j}(t)$ in (\arabicsection.\arabicequation). Now we will work with $\bm{\phi}_{j}(t).$

According to (Chen, 2007, Section 2.3), we have that for some pre-chosen orthonormal basis functions on $[0,1]$ , denoted as $\{\alpha_{k}(t)\}$

\bm{\phi}_{j,l}(t):=\sum_{k=1}^{c}a_{jk,l}\alpha_{k}(t)+\mathrm{O}(c^{-d}),\ 1% \leq l\leq j,

(\arabicsection.\arabicequation)

where $c$ is the number of basis functions. Here $\{a_{jk,l}\}$ are the coefficients to be estimated. In fact, using (Ding & Zhou, 2023, Theorem 2.11) and a discussion similar to (3.8) therein, we have by (\arabicsection.\arabicequation) that

x_{i}=\sum_{l=1}^{j}\sum_{k=1}^{c}a_{jk,l}z_{kl}+\epsilon_{i,j}+\mathrm{O}_{% \ell^{2}}\left(j^{2}/n+jc^{-d}\right),\ i>j,

(\arabicsection.\arabicequation)

where $\epsilon_{i,j}$ is defined in (\arabicsection.\arabicequation) and $z_{kl}(i/n)=\alpha_{k}(i/n)x_{i-l}$ . Observe that a key component in (\arabicsection.\arabicequation) is the approximation of $\rho_{j}(\cdot)$ by local best linear forecast coefficients as in (\arabicsection.\arabicequation).

Using (\arabicsection.\arabicequation), we can estimate all the coefficients $a_{jk,l}$ ’s using one OLS regression. In particular, we stack $a_{jk,l},1\leq l\leq j,1\leq k\leq c$ as a vector $\bm{\beta}_{j}\in\mathbb{R}^{jc},$ then the OLS estimator for $\bm{\beta}_{j}$ can be written as $\widehat{\bm{\beta}}_{j}:=(Y_{j}^{*}Y_{j})^{-1}Y_{j}^{*}\bm{x}_{j},$ where $Y_{j}^{*}\in\mathbb{R}^{jc\times(n-j)}$ is the design matrix of (\arabicsection.\arabicequation) and $\bm{x}_{j}=(x_{j+1},\cdots,x_{n})^{*}\in\mathbb{R}^{n-j}.$

After estimating the $a_{jk,l}$ ’s, $\rho_{j}(t)\equiv\bm{\phi}_{j,j}(t)$ is estimated using (\arabicsection.\arabicequation) as

\widehat{\rho}_{j}(t)=\widehat{\bm{\beta}}_{j}^{*}\mathbb{B}_{j,j}(t),

(\arabicsection.\arabicequation)

where $\mathbb{B}_{j,l}(t)\in\mathbb{R}^{jc}$ has $j$ blocks and the $l$ th block is $\mathbf{B}(t):=(\alpha_{1}(t),\cdots,\alpha_{c}(t))^{*}\in\mathbb{R}^{c},1\leq l\leq j$ and zeros otherwise.

\arabicsection.\arabicsubsection Multiplier bootstrap based inference

In this subsection we propose a multiplier bootstrap procedure to infer the PACFs. Statistical inference of the PACFs plays an important role in stationary time series analysis. For example, it can be used to determine the order of an AR process and check whether the time series (or residuals after an ARIMA model fitting) is white noise. We refer the readers to Chapter 3 of Shumway & Stoffer (2017) for more details. However, the analogs for locally stationary time series are largely missing. We aim to fill the gap in this section.

Based on our estimators in (\arabicsection.\arabicequation), we can conduct various important tests on $\rho_{j}(t)$ in (\arabicsection.\arabicequation). For example, we can test whether the PACFs are identical to functions of interest that $\rho_{j}(t)=f_{j}(t)$ for some given functions $f_{j}(t).$ For another instance, we can check whether a group of the PACFs are time-invariant that $\rho_{j}(t)\equiv\rho_{j},\ h_{1}\leq j\leq h_{2}$ for some integers $h_{1}$ and $h_{2},$ where $\rho_{j}(t)\equiv\rho_{j}$ means $\rho_{j}(t)=\rho_{j}$ for all $t$ . Note that if we set $\rho_{j}=0,$ it reduces to the testing the significance of the PACFs. Especially, when $h_{1}=1$ and $h_{2}=\infty,$ it is equivalent to test whether the time series is white noise. While we are able to conduct several different tests on the PACFs, in this paper, motivated by their applications in model selection and goodness of fit, for conciseness, we will focus on two important such tests.

First, we are interested in testing

\mathbf{H}_{01}:\ \rho_{j}(t)\equiv 0\ \ \text{vs}\ \ \mathbf{H}_{a1}:\ \rho_{% j}(t)\not\equiv 0,\ \text{for some}\ j\geq 1.

(\arabicsection.\arabicequation)

The hypothesis (\arabicsection.\arabicequation) tests the significance of a single PACF. Similar to the stationary setting in Brockwell & Davis (1987), it can be used to select the order of a locally stationary AR process; see Remark B.\arabicthm of our supplement for more discussions.

Second, we are also interested in testing the significance for all the lags that

\mathbf{H}_{02}:\rho_{j}(t)\equiv 0\ \text{for all}\ j\geq 1\ \ \text{vs}\ \ % \mathbf{H}_{a2}:\rho_{k}(t)\not\equiv 0\ \text{for some}\ k\geq 1.

(\arabicsection.\arabicequation)

The hypothesis (\arabicsection.\arabicequation) tests white noise (or lack of serial correlation) of the underlying locally time series. For stationary white noise, the well-known Box–Pierce (BP) test statistic Box & Pierce (1970) with fixed lag truncation number is probably the most commonly used statistic. Later on, such a test was extended to locally stationary white noise in Goerg (2012). We emphasize that the Portmanteau-type BP tests in Box & Pierce (1970); Goerg (2012) used the autocorrelation functions (ACFs) instead of the PACFs. However, the estimation of ACFs for locally stationary time series involves the estimation of the time-varying marginal variances which requires the choice of additional tuning parameters and may lead to deteriorated estimation accuracy in finite samples. Inspired by the above challenges and the discussions of Section 9.4 of Brockwell & Davis (1987), we will propose a PACF-based Portmanteau test.

We mention again that as discussed in Remark \arabicsection.\arabicthm, when $\tau$ is large, we have that $\sup_{t}\rho_{j}(t)=\mathrm{o}(n^{-1/2})$ for $j\gg j^{*}$ defined in (\arabicsection.\arabicequation). Therefore, from an inferential viewpoint, $\rho_{j}(\cdot)$ for $j\gg j^{*}$ can be effectively treated as zero. Consequently, we only need to consider the setting $j\leq j^{*}$ for (\arabicsection.\arabicequation) and (\arabicsection.\arabicequation) once $j^{*}$ is identified.

\arabicsection.\arabicsubsection.\arabicsubsubsection Test statistics

In this section, we propose the test statistics for (\arabicsection.\arabicequation) and (\arabicsection.\arabicequation). First, when the null hypothesis in (\arabicsection.\arabicequation) holds, the following statistic $T_{1}\equiv T_{1}(j)$ should be small

T_{1}\equiv T_{1}(j):=\int_{0}^{1}\widehat{\rho}_{j}^{2}(t)\mathrm{d}t.

(\arabicsection.\arabicequation)

Therefore, it is natural to use (\arabicsection.\arabicequation) to test (\arabicsection.\arabicequation).

Second, to test (\arabicsection.\arabicequation), motivated by the BP test Box & Pierce (1970), we may want to directly use the following statistic

T_{\text{BP}}\equiv T_{\text{BP}}(\mathsf{h}):=\sum_{k=1}^{\mathsf{h}}\int_{0}% ^{1}\widehat{\rho}_{k}(t)^{2}\mathrm{d}t,\ \text{for some large}\ \mathsf{h}% \geq j^{*},

(\arabicsection.\arabicequation)

where we recall $j^{*}$ in (\arabicsection.\arabicequation). Even though it is natural to use $T_{\text{BP}}$ , as described in Section \arabicsection.\arabicsubsection, in order to obtain its value, we need to do $\mathsf{h}$ high dimensional OLS regressions which can be computationally expensively, especially when $\mathsf{h}$ is large (or equivalently, $\tau$ is small). To address this issue, we consider an $\mathsf{h}$ order best linear prediction as in (\arabicsection.\arabicequation). That is, $x_{i}=\sum_{k=1}^{\mathsf{h}}\phi_{ik,\mathsf{h}}x_{i-k}+\epsilon_{i,\mathsf{h% }}.$ According to Theorem 2.11 of Ding & Zhou (2023), by setting $j=\mathsf{h}$ in (\arabicsection.\arabicequation), we find that

x_{i}=\sum_{k=1}^{\mathsf{h}}\phi_{\mathsf{h},k}(i/n)x_{i-k}+\epsilon_{i,% \mathsf{h}}+\mathrm{O}_{\ell^{2}}\left(\mathsf{h}^{2}/n\right).

(\arabicsection.\arabicequation)

The smooth coefficients $\{\phi_{\mathsf{h},k}(\cdot)\}$ can be estimated via the sieve method using only one high-dimensional OLS as in Section \arabicsection.\arabicsubsection whose estimators are denoted as $\{\widehat{\phi}_{\mathsf{h},k}(\cdot)\}.$

Now we define another statistic

T_{2}\equiv T_{2}(\mathsf{h}):=\sum_{k=1}^{\mathsf{h}}\int_{0}^{1}\widehat{% \phi}_{\mathsf{h},k}^{2}(t)\mathrm{d}t.

(\arabicsection.\arabicequation)

As will be seen later in Theorem \arabicsection.\arabicthm below, when (\arabicsection.\arabicequation) holds, under some mild conditions on $\mathsf{h}$ (c.f. (\arabicsection.\arabicequation)), $T_{2}$ will be close to $T_{\text{BP}}$ while the calculation of $T_{2}$ only needs one OLS regression. Therefore, we will use (\arabicsection.\arabicequation) to test (\arabicsection.\arabicequation).

\arabicsection.\arabicsubsection.\arabicsubsubsection Practical implementation

We point out that it is still difficult to directly use $T_{1}$ and $T_{2}$ since the variances in their limiting Gaussian distributions are usually hard to estimate and hence plug-in estimators are unavailable; see Theorems \arabicsection.\arabicthm and \arabicsection.\arabicthm for more details. To address this issue, we utilize the multiplier bootstrap procedure as in Ding & Zhou (2023); Zhou (2013a). We first explain how to construct the bootstrapped statistics. Using the sieve estimates as in Section \arabicsection.\arabicsubsection, denote the residual

\widehat{\epsilon}_{i,\ell}:=x_{i}-\sum_{k=1}^{\ell}\widehat{\phi}_{k}(i/n)x_{% i-k},

(\arabicsection.\arabicequation)

where $\ell=j$ for $T_{1}$ and $\ell=\mathsf{h}$ for $T_{2}.$

Let $\widehat{\bm{w}}_{i}=\bm{x}_{\ell,i}\widehat{\epsilon}_{i,\ell},$ where $\bm{x}_{\ell,i}=(x_{i-1},\cdots,x_{i-\ell})^{*}.$ Given a block size $m,$ we denote $\widehat{\Phi}\equiv\widehat{\Phi}(\ell,m)\in\mathbb{R}^{\ell c}$ as

\widehat{\Phi}:=\frac{1}{\sqrt{n-m-\ell+1}\sqrt{m}}\sum_{i=\ell+1}^{n-m}\left[% \left(\sum_{j=i}^{i+m}\widehat{\bm{w}}_{\ell,j}\right)\otimes\left(\mathbf{B}(% i/n)\right)\right]R_{i},

(\arabicsection.\arabicequation)

where $\otimes$ is the Kronecker product and $R_{i},\ell+1\leq i\leq n-m$ are i.i.d. standard Gaussian random variables which are independent of the observed time series. Recall the discussions around (\arabicsection.\arabicequation) and $Y_{\ell}^{*}$ is the design matrix. Denote $\widehat{\Sigma}_{1}:=\frac{1}{n}Y_{j}^{*}Y_{j}$ and $\widehat{\Sigma}_{2}:=\frac{1}{n}Y_{\mathsf{h}}^{*}Y_{\mathsf{h}}.$ Moreover, let $\mathbf{M}\in\mathbb{R}^{\ell c\times\ell c}$ be a diagonal block matrix whose only non-zero part is the identity matrix lies in the last diagonal block. Inspired by Remark \arabicsection.\arabicthm below, we use the following statistics $\widehat{\mathcal{T}}_{k},k=1,2,$ to mimic the distribution of $nT_{k},k=1,2.$ Let $\widehat{\Phi}_{1}$ be constructed as in (\arabicsection.\arabicequation) using $\ell=j$ and $\widehat{\Phi}_{2}$ be constructed using $\ell=\mathsf{h}.$ Then we denote

\widehat{\mathcal{T}}_{1}=\widehat{\Phi}_{1}^{*}\widehat{\Sigma}^{-1}\mathbf{M% }\widehat{\Sigma}_{1}^{-1}\widehat{\Phi}_{1},\ \widehat{\mathcal{T}}_{2}=% \widehat{\Phi}_{2}^{*}\widehat{\Sigma}_{2}^{-2}\widehat{\Phi}_{2}.

(\arabicsection.\arabicequation)

We point out that for the implementation, one needs to select some large value of $\mathsf{h}$ to construct $T_{2}.$ We discuss how to choose this parameter in Section C of our supplement.

Finally, based on the above result, we propose the following Algorithm 1 for the practical implementation. Note that in order to implement Algorithm 1, two tuning parameters, the number of basis functions $c$ and the block size $m$ , have to been chosen properly. In our $\mathtt{R}$ package $\mathtt{Sie2nts}$ , these parameters can be chosen automatically using the function $\mathtt{auto.pacf.test}$ according to the methods provided in Section C of our supplement Ding & Zhou (2024).

Algorithm 1 Multiplier Bootstrap

Inputs: The lag $j$ for (\arabicsection.\arabicequation) or $\mathsf{h}$ for (\arabicsection.\arabicequation), type I error rate $\alpha,$ tuning parameters $c$ and $m$ chosen by the data-driven procedure demonstrated in Section C of our supplement, time series $\{x_{i}\},$ and sieve basis functions.

Step one: Compute $\widehat{\Sigma}_{1}^{-1}$ using $n(Y_{j}^{*}Y_{j})^{-1}$ for (\arabicsection.\arabicequation) or $\widehat{\Sigma}_{2}^{-1}$ using $n(Y_{\mathsf{h}}^{*}Y_{\mathsf{h}})^{-1}$ for (\arabicsection.\arabicequation), and the residuals $\{\widehat{\epsilon}_{i,\ell}\}_{i=\ell+1}^{n}$ according to (\arabicsection.\arabicequation) with $\ell=j$ for (\arabicsection.\arabicequation) or $\ell=\mathsf{h}$ for (\arabicsection.\arabicequation).

Step two: Generate $B$ (say 1,000) i.i.d. copies of $\{\widehat{\Phi}^{(s)}\}_{s=1}^{B}$ according to (\arabicsection.\arabicequation). Compute $\widehat{\mathcal{T}}_{k}^{s}$ , $k=1,2$ , $s=1,2,\cdots,B,$ correspondingly as in (\arabicsection.\arabicequation).

Step three: Let $\widehat{\mathcal{T}}^{(1)}_{k}\leq\widehat{\mathcal{T}}^{(2)}_{k}\leq\cdots% \leq\widehat{\mathcal{T}}^{(B)}_{k}$ be the order statistics of $\widehat{\mathcal{T}}^{s}_{k},s=1,2,\cdots,B.$ Reject $\mathbf{H}_{01}$ in (\arabicsection.\arabicequation) at the level $\alpha$ if $nT_{1}>\widehat{\mathcal{T}}_{1}^{(\lfloor B(1-\alpha)\rfloor)},$ where $\lfloor x\rfloor$ denotes the largest integer smaller or equal to $x.$ Reject $\mathbf{H}_{02}$ in (\arabicsection.\arabicequation) at the level $\alpha$ if $nT_{2}>\widehat{\mathcal{T}}_{2}^{(\lfloor B(1-\alpha)\rfloor)}.$ Let $B_{k}^{*}=\max\{r:\widehat{\mathcal{T}}^{r}_{k}\leq nT_{k}\},k=1,2.$

Output: $p$ -value of the tests (\arabicsection.\arabicequation) and (\arabicsection.\arabicequation) can be computed, respectively as $1-\frac{B^{*}_{k}}{B},k=1,2.$

\arabicsection Theoretical analysis

In this section, we provide some theoretical analysis on our estimation and inference procedures. Till the end of the paper, for notational simplicity, we assume that the locally stationary time series admits the general physical representation equipped with the physical dependence measures (see (B.\arabicequation) and (B.\arabicequation) of our supplement). In addition, we need the following assumptions.

Assumption \arabicsection.\arabicthm.

Throughout the paper, we suppose the followings holds:

(1).

For all sufficiently large $n\in\mathbb{N},$ we assume that there exists a universal constant $\kappa>0$ that

\lambda_{n}(\operatorname{Cov}(x_{1},\cdots,x_{n}))\geq\kappa,

(\arabicsection.\arabicequation)

where $\lambda_{n}(\cdot)$ is the smallest eigenvalue of the given matrix and $\operatorname{Cov}(\cdot)$ is the covariance matrix of the given vector.

(2).

For all $n\in\mathbb{N},$ $1\leq k\leq n$ and $-k+1\leq r\leq n-k,$ we assume that there exists some constant $\tau>1$ such that (\arabicsection.\arabicequation) holds. In addition, we assume that $\sup_{i,n}\mathbb{E}|x_{i}|<\infty.$

(3).

For some given integer $d>0$ , we assume that for any $j\geq 0$

\gamma(t,j)\in C^{d}([0,1]).

(\arabicsection.\arabicequation)

The conditions in the above assumption are mild and can be satisfied by many commonly used locally stationary time series. Due to space constraint, we leave some discussions to Remark B.\arabicthm of supplement Ding & Zhou (2024).

\arabicsection.\arabicsubsection Uniform consistency

In what follows, we establish the consistency for our estimators. Recall $\mathbf{B}(t)$ below (\arabicsection.\arabicequation). Denote

\xi_{c}=\sup_{t}\sup_{1\leq i\leq c}|\alpha_{i}(t)|,\ \zeta_{c}=\sup_{t}|% \mathbf{B}(t)|.

(\arabicsection.\arabicequation)

The following mild assumption will be needed to ensure a consistent estimation, which has been used frequently in the literature, see Ding & Zhou (2020, 2023, 2021); Vogt (2012). Recall $\gamma(\cdot,\cdot)$ in (\arabicsection.\arabicequation). For all $j=1,2,\cdots,j_{*},$ denote $\Sigma^{j}(t)\in\mathbb{R}^{j\times j}$ whose $(k,l)$ th entry is $\Sigma_{kl}^{(j)}(t)=\gamma(t,|k-l|).$

Assumption \arabicsection.\arabicthm.

For $j=1,2,3,\cdots,j_{*},$ denote the long-run integrated covariance matrix as

\Sigma^{(j)}=\int_{0}^{1}\Sigma^{(j)}(t)\otimes\left(\mathbf{B}(t)\mathbf{B}(t% )^{*}\right)\mathrm{d}t,

(\arabicsection.\arabicequation)

where we recall that $\otimes$ is the Kronecker product. We assume that the eigenvalues of $\Sigma^{(j)}$ are bounded above and also away from zero by some universal constants.

Then we proceed to state the main results of this section. Recall (\arabicsection.\arabicequation). Denote

\Psi(j,c):=j\xi_{c}\zeta_{c}\sqrt{\frac{c}{n}}\left(1+\frac{j^{2}}{n}+jc^{-d}% \right)+c^{-d}.

(\arabicsection.\arabicequation)

Theorem \arabicsection.\arabicthm.

Suppose Assumptions \arabicsection.\arabicthm and \arabicsection.\arabicthm hold true. Moreover, for $j\leq j^{*}$ satisfying

jc\left(\frac{\xi_{c}^{2}}{\sqrt{n}}+\frac{\xi_{c}^{2}n^{\frac{2}{\tau+1}}}{n}% \right)=\mathrm{o}(1),

(\arabicsection.\arabicequation)

we have that for our estimator (\arabicsection.\arabicequation)

\sup_{t}\left|\widehat{\rho}_{j}(t)-\rho_{j}(t)\right|=\mathrm{O}_{\mathbb{P}}% \left(\Psi(j,c)\right).

(\arabicsection.\arabicequation)

Remark \arabicsection.\arabicthm.

Theorem \arabicsection.\arabicthm implies that our proposed estimator (\arabicsection.\arabicequation) is uniformly consistent under mild conditions. First, the condition (\arabicsection.\arabicequation) ensures that $n^{-1}Y_{j}^{*}Y_{j}$ in the OLS estimator will convergence to $\Sigma^{(j)}$ which guarantees the regular behavior of $\widehat{\bm{\beta}}_{j}.$ In fact, (\arabicsection.\arabicequation) can be easily satisfied. Note that $\xi_{c}$ and $\zeta_{c}$ can be calculated for the specific sieve basis functions and so does the convergence rate in (\arabicsection.\arabicequation). For example, when $\{\alpha_{k}(t)\}$ are Fourier basis functions or normalized orthogonal polynomials, $\xi_{c}=\mathrm{O}(1)$ and $\zeta_{c}=\mathrm{O}(\sqrt{c}).$ Consequently, when $\xi_{c}=\mathrm{O}(1),$ even for $j=j^{*}$ in (\arabicsection.\arabicequation), (\arabicsection.\arabicequation) only requires that $c(n^{-1/2+1/(2(\tau-1))+n^{-1+1/(2(\tau-1))+2/(\tau+1)}})\ll 1.$ In particular, if we set $c=\mathrm{O}(n^{a})$ for some small constant $0<a<1/2,$ we only need $\tau>1+\frac{1}{1-2a}.$

Second, for the rate $\Psi(j,c)$ in (\arabicsection.\arabicequation) and (\arabicsection.\arabicequation), when $\xi_{c}=\mathrm{O}(1)$ and $\zeta_{c}=\mathrm{O}(\sqrt{c}),$ even for $j=j^{*},$ it reads $c^{-d}(1+cn^{-1/2+1/(\tau-1)})+cn^{-1/2+1/(2(\tau-1))}+cn^{-3/2+3/(2(\tau-1))}$ . Therefore, for sufficiently large $d$ and $\tau,$ it has an order of $cn^{-1/2+\epsilon},$ for some small constant $\epsilon>0.$

\arabicsection.\arabicsubsection Asymptotic normality and power analysis for the proposed statistics

In this section, we study the accuracy and power of the proposed statistics $T_{1}$ in (\arabicsection.\arabicequation) and $T_{2}$ in (\arabicsection.\arabicequation). We first prepare some notations. Recall (\arabicsection.\arabicequation). Following the conventions below (\arabicsection.\arabicequation), for $\ell=j$ regarding $T_{1}$ and $\ell=\mathsf{h}$ regarding $T_{2},$ we denote

\bm{w}_{i}=\bm{x}_{\ell,i}\epsilon_{i,\ell},\ i>\ell,

(\arabicsection.\arabicequation)

where we recall $\bm{x}_{\ell,i}=(x_{i-1},\cdots,x_{i-\ell})^{*}\in\mathbb{R}^{\ell}.$ According to (3.17) of Ding & Zhou (2023) or Lemma 3.1 of Ding & Zhou (2021), we see that $\bm{w}_{i}$ has a physical representation in the sense that for some measurable function $\mathbf{V},$ we have

\bm{w}_{i}=\mathbf{V}(\frac{i}{n},\mathcal{F}_{i}),\ i>\ell,

(\arabicsection.\arabicequation)

where $\mathcal{F}_{i}=(\cdots,\eta_{i-1},\eta_{i})$ and $\eta_{i},\ i\in\mathbb{Z}$ are i.i.d centered random variables. Denote the long-run covariance matrix of $\bm{w}_{i}$ as $\Pi(t)=\sum_{j=-\infty}^{\infty}\text{Cov}\Big{(}\mathbf{V}(t,\mathcal{F}_{0})% ,\mathbf{V}(t,\mathcal{F}_{j})\Big{)},$ and the integrated long-run covariance matrix as $\Pi=\int_{0}^{1}\Pi(t)\otimes(\mathbf{B}(t)\mathbf{B}^{*}(t))\mathrm{d}t.$ In what follows, to ease our discussion, we assume that $c$ is of the form

c=\mathrm{O}(n^{a}),\ 0<a<1.

(\arabicsection.\arabicequation)

Armed with the above notations, we now proceed to provide the theoretical properties of the statistic $T_{1}.$ Recall the matrix $\mathbf{M}$ in (\arabicsection.\arabicequation). For $s\in\mathbb{N},$ define

g_{1,s}=(\text{Tr}[\Pi^{1/2}\overline{\Sigma}_{01}^{-1}\mathbf{M}\overline{% \Sigma}_{01}^{-1}\Pi^{1/2}]^{s})^{1/s},\ \overline{\Sigma}_{01}=\begin{pmatrix% }\mathbf{I}_{c}&\bm{0}\\ \bm{0}&\Sigma_{01}\end{pmatrix},

(\arabicsection.\arabicequation)

where $\Sigma_{01}:=\int_{0}^{1}\Sigma^{j}(t)\otimes(\mathbf{B}(t)\mathbf{B}^{*}(t))% \mathrm{d}t$ and $\Pi$ is defined from (\arabicsection.\arabicequation) using $\ell=j.$

Theorem \arabicsection.\arabicthm.

Suppose the assumptions of Theorem \arabicsection.\arabicthm hold. Then we have that

\arabicenumi.

Suppose Assumption A.\arabicthm of our supplement holds and

\Psi(j,c)=\mathrm{o}(1),

(\arabicsection.\arabicequation)

we have that when $\mathbf{H}_{01}$ in (\arabicsection.\arabicequation) holds

\mathbb{T}_{1}:=\frac{nT_{1}-g_{1,1}}{g_{1,2}}\Rightarrow\mathcal{N}(0,2).

(\arabicsection.\arabicequation)

\arabicenumi.

When $\mathbf{H}_{a1}$ in (\arabicsection.\arabicequation) holds in the sense that

\int_{0}^{1}\rho_{j}(t)^{2}\mathrm{d}t>\mathfrak{C},\ \text{where}\ \mathfrak{% C}:=C_{\alpha}\frac{\sqrt{jc}}{n},

(\arabicsection.\arabicequation)

where $C_{\alpha}\equiv C_{\alpha}(n)\rightarrow\infty$ as $n\rightarrow\infty,$ assuming (\arabicsection.\arabicequation), then we have that for any $\alpha\in(0,1)$

\mathbb{P}\left(\left|\mathbb{T}_{1}\right|\geq\sqrt{2}\mathcal{Z}_{1-\alpha}% \right)\rightarrow 1,\ n\rightarrow\infty,

where $\mathcal{Z}_{1-\alpha}$ is the $(1-\alpha)$ th quantile of the standard Gaussian distribution.

The above theorem establishes the asymptotic normality for our proposed statistic $T_{1}$ in (\arabicsection.\arabicequation) concerning (\arabicsection.\arabicequation). Moreover, (\arabicsection.\arabicequation) shows that our proposed statistic can have asymptotic power one under weak local alternative. Then we provide the theoretical properties of the statistic $T_{2}.$ For $s\in\mathbb{N},$ define

g_{2,s}=(\text{Tr}[\Pi^{1/2}\overline{\Sigma}_{02}^{-2}\Pi^{1/2}]^{s})^{1/s},% \ \overline{\Sigma}_{02}=\begin{pmatrix}\mathbf{I}_{c}&\bm{0}\\ \bm{0}&\Sigma_{02}\end{pmatrix},

(\arabicsection.\arabicequation)

where $\Sigma_{02}:=\int_{0}^{1}\Sigma^{\mathsf{h}}(t)\otimes(\mathbf{B}(t)\mathbf{B}% ^{*}(t))\mathrm{d}t$ and $\Pi$ is defined from (\arabicsection.\arabicequation) using $\ell=\mathsf{h}.$

Theorem \arabicsection.\arabicthm.

Suppose the assumptions of Theorem \arabicsection.\arabicthm hold. Then we have that

\arabicenumi.

Suppose Assumption A.\arabicthm of our supplement holds and

\mathsf{h}\left(\frac{\mathsf{h}^{4}}{n}+\mathsf{h}^{2}\Psi^{2}(\mathsf{h},c)% \right)=\mathrm{o}(1),

(\arabicsection.\arabicequation)

when $\mathbf{H}_{02}$ in (\arabicsection.\arabicequation) holds, we have that

\mathbb{T}_{2}:=\frac{nT_{2}-g_{2,1}}{g_{2,2}}\Rightarrow\mathcal{N}(0,2).

(\arabicsection.\arabicequation)

\arabicenumi.

Recall (\arabicsection.\arabicequation). Suppose (\arabicsection.\arabicequation) holds. Then when $\mathbf{H}_{02}$ in (\arabicsection.\arabicequation) holds, we have that

T_{\text{BP}}=T_{2}+\mathrm{O}_{\mathbb{P}}\left(\mathsf{h}\Psi^{2}(\mathsf{h}% ,c)+\frac{\mathsf{h}^{5}}{n^{2}}+\sum_{k=1}^{h_{2}}\Psi^{2}(k,c)\right).

(\arabicsection.\arabicequation)

Consequently, if we further assume that

\frac{n\left[\mathsf{h}\Psi^{2}(\mathsf{h},c)+\frac{\mathsf{h}^{5}}{n^{2}}+% \sum_{k=1}^{\mathsf{h}}\Psi^{2}(k,c)\right]}{\sqrt{\mathsf{h}c}}=\mathrm{o}(1),

(\arabicsection.\arabicequation)

then

\frac{nT_{\text{BP}}-g_{2,1}}{g_{2,2}}\Rightarrow\mathcal{N}(0,2).

(\arabicsection.\arabicequation)

\arabicenumi.

When $\mathbf{H}_{a2}$ in (\arabicsection.\arabicequation) holds in the sense that

\sum_{k=1}^{\mathsf{h}}\int_{0}^{1}\rho_{k}^{2}(t)\mathrm{d}t>\mathfrak{C},\ % \text{where}\ \mathfrak{C}:=C_{\alpha}\frac{\sqrt{\mathsf{h}c}}{n},

(\arabicsection.\arabicequation)

where $C_{\alpha}\equiv C_{\alpha}(n)\rightarrow\infty$ as $n\rightarrow\infty,$ assuming the assumptions of parts 1 and 2 hold, then we have that for any $\alpha\in(0,1)$

\mathbb{P}\left(\left|\mathbb{T}_{2}\right|\geq\sqrt{2}\mathcal{Z}_{1-\alpha}% \right)\rightarrow 1,\ n\rightarrow\infty,

where $\mathcal{Z}_{1-\alpha}$ is the $(1-\alpha)$ th quantile of the standard Gaussian distribution.

Theorem \arabicsection.\arabicthm implies that $T_{\text{BP}}$ and $T_{2}$ have the same Gaussian distribution and power performance asymptotically so that we can directly use $T_{2}$ which only requires one OLS regression for estimation. We point out that (\arabicsection.\arabicequation) and (\arabicsection.\arabicequation) basically impose some upper bound conditions for $\mathsf{h}.$ As discussed in Remark \arabicsection.\arabicthm, if one chooses Fourier or orthogonal polynomial as the basis functions, when $d$ and $\tau$ are large enough, we require $\mathsf{h}\ll n^{1/5}$ to guarantee (\arabicsection.\arabicequation). Analogously, (\arabicsection.\arabicequation) requires that $\mathsf{h}\ll\min\{c^{1/5},n^{1/5}\}.$ Recall (\arabicsection.\arabicequation). In particular, $j^{*}=o(n^{1/5})$ when $\tau>2.5$ . Note that $\mathsf{h}$ is only required to be larger or equal to $j^{*}$ . As a result, the above constraints on $\mathsf{h}$ is not restrictive.

Remark \arabicsection.\arabicthm.

We add a remark on the assumptions of the parameters. First, (\arabicsection.\arabicequation) and (\arabicsection.\arabicequation) are mainly used to guarantee that $nT_{k},k=1,2,$ can be written as a quadratic form in terms of $\bm{w}_{i}$ in (\arabicsection.\arabicequation) under the null hypotheses. More specifically, denote $\mathbf{X}=\frac{1}{\sqrt{n}}\sum_{i=\ell+1}^{n}(\bm{w}_{i}\otimes\mathbf{B}(i% /n)),$ we can conclude from Ding & Zhou (2023) that

nT_{1}=\mathbf{X}^{*}\overline{\Sigma}_{01}^{-1}\mathbf{M}\overline{\Sigma}_{0% 1}^{-1}\mathbf{X}+\mathrm{o}_{\mathbb{P}}(1),\ \ nT_{2}=\mathbf{X}^{*}% \overline{\Sigma}_{02}^{-2}\mathbf{X}+\mathrm{o}_{\mathbb{P}}(1).

(\arabicsection.\arabicequation)

Second, (\arabicsection.\arabicequation) guarantees the difference between $T_{\text{BP}}$ and $T_{2}$ is negligible when $\mathbf{H}_{02}$ in (\arabicsection.\arabicequation) holds.

Before concluding this section, we summarize the properties of $\widehat{\mathcal{T}}_{k},k=1,2,$ in (\arabicsection.\arabicequation) which explains the validity for the bootstrap procedure. The motivation comes from the arguments as in Remark \arabicsection.\arabicthm that the statistic $nT_{k},k=1,2,$ are essentially quadratic forms of the locally stationary vector (c.f. (\arabicsection.\arabicequation)) and the statistics are asymptotically Gaussian. Consequently, it is possible to mimic their asymptotic distribution using a multiplier bootstrap procedure.

Corollary \arabicsection.\arabicthm.

Suppose the assumptions of Theorems \arabicsection.\arabicthm and \arabicsection.\arabicthm hold. Moreover, assume the time series has finite fourth moment and $m\rightarrow\infty$ as $n\rightarrow\infty$ and

\sqrt{\mathsf{h}}\zeta_{c}^{2}c^{-1/2}\left(\frac{m}{n}+\frac{1}{m}\right)=% \mathrm{o}(1).

Then there exists a sequence of sets $\mathcal{A}_{n}$ such that $\mathbb{P}(\mathcal{A}_{n})=1-\mathrm{o}(1)$ and under the event $\mathcal{A}_{n},$ we have that conditional on the data $\{x_{i}\},$ the results in Theorems \arabicsection.\arabicthm and \arabicsection.\arabicthm still hold by replacing $nT_{k}$ with $\widehat{\mathcal{T}}_{k},k=1,2.$

\arabicsection Numerical simulations and real data analysis

In this section, we use some numerical simulations and a real data analysis to illustrate the usefulness of our estimation and inference procedures of the PACFs. All the calculations, implementations and plots can be done using a few lines of coding with our $\mathtt{R}$ package $\mathtt{Sie2nts}.$

\arabicsection.\arabicsubsection Numerical simulations

In this section, we conduct numerical simulations to illustrate the usefulness of our methodologies using both stationary and non-stationary models. Due to space constraint, we focus on reporting the results of AR type models in the following. Additional simulation results on other types of models can be found in Section D of our supplement.

In what follows, for some $\delta_{1},\delta_{2},\in[0,0.5],$ we consider the stationary AR(2) process

x_{i}=\delta_{1}x_{i-1}+\delta_{2}x_{i-2}+\epsilon_{i},

(\arabicsection.\arabicequation)

and the locally stationary AR(2) process

x_{i}=\delta_{1}\sin(2\pi i/n)x_{i-1}+\delta_{2}\cos(2\pi i/n)x_{i-2}+\left(0.% 4+0.4|\sin(2\pi i/n)|\right)\epsilon_{i},

(\arabicsection.\arabicequation)

where $\epsilon_{i},1\leq i\leq n,$ are i.i.d. standard Gaussian random variables. Note that when $\delta_{1}=\delta_{2}=0,$ (\arabicsection.\arabicequation) reduces to a standard white noise process and (\arabicsection.\arabicequation) reduces to a time-varying white noise process.

\arabicsection.\arabicsubsection.\arabicsubsubsection Estimation of PACFs

We first estimate the PACFs using the sieve method as introduced in Section \arabicsection.\arabicsubsection. For concreteness and due to space constraint, regarding $\delta_{1}$ and $\delta_{2}$ in (\arabicsection.\arabicequation) and (\arabicsection.\arabicequation), we only report the results for the choice $\delta_{1}=0.5$ and $\delta_{2}=0.3.$ Note that similar results and conclusions can also be made for other choices.

In Figure \arabicfigure, we estimate the PACFs of the first 10 lags for both models (\arabicsection.\arabicequation) and (\arabicsection.\arabicequation) using the estimators in (\arabicsection.\arabicequation). We use the Legendre polynomials as the basis functions and the number $c$ can be chosen using the cross validation method as described in Section C. The computations of the PACFs can be obtained directly using the function $\mathtt{sie.plot.pacf}$ from our $\mathtt{R}$ package. From these plots, we can see that our method applies to both models and obtain reasonably accurate estimates. According to the cut-off properties of the PACFs of the AR models, these plots have also implied that the time series may be generated from some AR(2) models. In Section D.\arabicsubsection of our supplement, we compare our proposed method with the ones introduced in Killick et al. (2020) and find that our method is generally more accurate than Killick et al. (2020) in terms of mean integrated squared error. The main reason is that the sieve method is adaptive to the smoothness of the covariance structure and has less boundary issues compared with the kernel method. In addition, more simulation results on other types of models can be found in Section D.\arabicsubsection of our supplement and similar conclusions can also be made.

Refer to caption — (a) PACFs for model (\arabicsection.\arabicequation).

\arabicsection.\arabicsubsection.\arabicsubsubsection Inference of PACFs

In this section, we examine the accuracy and sensitivity of our proposed Algorithm 1 when applied to testing (\arabicsection.\arabicequation) and (\arabicsection.\arabicequation). We first investigate the accuracy. To test (\arabicsection.\arabicequation) for some individual PACF at lag $j$ , in the context of (\arabicsection.\arabicequation) and (\arabicsection.\arabicequation), we consider the following four settings: (1). $\delta_{1}=0.5,\delta_{2}=0,$ and $j=2$ ; (2). $\delta_{1}=0.5,\delta_{2}=0,$ and $j=4$ ; (3). $\delta_{1}=\delta_{2}=0.3,$ and $j=3;$ (4). $\delta_{1}=\delta_{2}=0.3,$ and $j=5.$ Moreover, to test (\arabicsection.\arabicequation) for the white noise, we consider the following setting in the context of (\arabicsection.\arabicequation) and (\arabicsection.\arabicequation): (5). $\delta_{1}=\delta_{2}=0.$

In Table \arabictable, we report the simulated type I error rates for all the above null settings for three different types of basis functions when $n=600.$ We can see that our Algorithm 1 is quite accurate for both tests (\arabicsection.\arabicequation) and (\arabicsection.\arabicequation).

	$\alpha=0.1$					$\alpha=0.05$
Basis/Setting	(1)	(2)	(3)	(4)	(5)	(1)	(2)	(3)	(4)	(5)
	Model (\arabicsection.\arabicequation)
Fourier	0.108	0.096	0.108	0.11	0.098	0.048	0.059	0.061	0.054	0.049
Legendre	0.109	0.1	0.096	0.103	0.108	0.06	0.053	0.055	0.064	0.048
Daubechies-9	0.096	0.102	0.11	0.101	0.095	0.058	0.061	0.054	0.058	0.043
	Model (\arabicsection.\arabicequation)
Fourier	0.098	0.11	0.108	0.099	0.107	0.048	0.062	0.057	0.048	0.039
Legendre	0.103	0.094	0.103	0.095	0.113	0.055	0.053	0.061	0.048	0.053
Daubechies-9	0.093	0.09	0.093	0.104	0.1	0.052	0.058	0.058	0.047	0.046

Table \arabictable: Simulated type I error rates. The results are reported based on 1,000 simulations. We can see that our multiplier bootstrap procedure is reasonably accurate for both

\alpha=0.1

and

\alpha=0.05

. In our

\mathtt{R}

package

\mathtt{Sie2nts},

we collect 31 different types basis functions. For visualization of various types of basis functions, we can use the

\mathtt{bs.gene}

and

\mathtt{bs.plot}

from our

\mathtt{R}

package

\mathtt{Sie2nts}.

The computations of the

p

-values can be obtained directly using the function

\mathtt{auto.pacf.test}

from the package.

Then we study the power. For definiteness, we report the results for the white noise test that

\mathbf{H}_{0}:\ \{x_{i}\}\ \text{is a white noise process}\ \ \text{Vs}\ \ % \mathbf{H}_{a}:\ \{x_{i}\}\ \text{is not a white noise process.}

(\arabicsection.\arabicequation)

More specifically, in view of the concerned models (\arabicsection.\arabicequation) and (\arabicsection.\arabicequation), $\mathbf{H}_{0}$ in (\arabicsection.\arabicequation) is equivalent to $\delta_{1}=\delta_{2}=0$ while $\mathbf{H}_{a}$ is an AR(1) alternative that $\delta_{2}=0,\delta_{1}>0.$ In Figure \arabicfigure below, we study the power of our proposed method when $\delta_{1}$ deviates away from zero. We can conclude that our proposed tests are reasonably powerful once the alternative deviates from the null.

Finally, for further visualization, using (\arabicsection.\arabicequation) as an example, in Figure \arabicfigure below, we provide two typical plots of the $p$ -values associated with each lag $j$ with (\arabicsection.\arabicequation) under both the null and alternative as in (\arabicsection.\arabicequation). Here for the alternative we set $\delta_{1}=0.5$ and $\delta_{2}=0$ in (\arabicsection.\arabicequation). From the plots we can easily distinguish the null and alternative. Moreover, it is clear that the plots suggest that the null is a white noise process and the alternative is an AR(1) process. Additionally, more simulation results on other types of models can be found in Section D.\arabicsubsection of our supplement and similar conclusions can also be made.

\arabicsection.\arabicsubsection Real data analysis

In this section, we apply our methods to the monthly Euro-Dollar exchange rate data set which has also been considered in Killick et al. (2020). The data set can be downloaded from EuroStat at https://ec.europa.eu/eurostat/web/products-datasets/-/ei_mfrt_m. We analyze the exchange rates from January 1999 until October 2017.

Following the traditions of financial data analysis, we consider the log returns of the exchange rate. Moreover, we apply our methods to estimate the PACFs and conduct inference on them. As can be seen from the left panel of Figure \arabicfigure, based on our inference, a white-noise-driven AR(1) model will be useful to model the data. This agrees with the findings in Killick et al. (2020). In fact, we can conduct the white noise test for the original time series as in (\arabicsection.\arabicequation) using Algorithm 1 and find the $p$ -value is $0.013.$ This shows that the time series is not a white noise. Furthermore, after fitting the time series with a time-varying AR(1) model (note that for AR(1) model, the coefficient is identical to $\widehat{\rho}_{1}(t)$ ), we can conduct the white noise test for the residuals and then find the $p$ -value is $0.54$ . This confirms that the AR(1) model is likely to be appropriate.

In the right panel of Figure \arabicfigure, we provide the estimation of the PACF of lag one which is also the coefficient of a time-varying AR(1) model. It can be seen that the even though the PACF is relatively stable, it experiences some smooth changes over time. In fact, one can directly generalize the significance test $T_{1}$ into a test of constancy by considering $T_{1}^{*}:=\int_{0}^{1}[\widehat{\rho}_{1}(t)-\int_{0}^{1}\widehat{\rho}_{1}(s% )\mathrm{d}s]^{2}\mathrm{d}t.$ The multiplier bootstrap of $T_{1}^{*}$ yields a $p$ -value of $0.031$ , indicating that $\rho_{1}(t)$ is likely to be time-varying. This concludes that a locally stationary AR(1) model may be more useful for the model fitting. Due to the pronounced responsiveness of Euro-Dollar exchange rates to the global economy, our plots serve as reasonably accurate reflections of prevailing global economic conditions. For instance, leading up to the global financial crisis (2005-2007), the PACF displays a consistent pattern. Subsequently, from 2008 onwards, there is a gradual decline observed, hitting its lowest point around 2013, coinciding with the recognized period of the financial crisis. Following this, the PACF demonstrates a resurgence. Such a visual representation offers valuable insights into the evolving dynamics of Euro-Dollar exchange rates, aiding in a deeper comprehension of their temporal behavior.

Supplementary file

In the supplement file, we provide the technical proofs, some additional remarks, practical methods for choosing the tuning parameters and additional simulation results.

A Technical proofs

A.\arabicsubsection Proofs of Section \arabicsection

In this subsection, we prove Theorem \arabicsection.\arabicthm. The strategies and ideas are similar to those of Theorems 2.4 and 2.11 of Ding & Zhou (2023). We focus on explaining the main differences.

{Proof}

[Proof of Theorem \arabicsection.\arabicthm] For the first part of the proof, it is analogous to the proof of equation (2.6) of Ding & Zhou (2023). We only sketch the proof. Due to similarity, we only prove the first control in (\arabicsection.\arabicequation). Recall (\arabicsection.\arabicequation). Set $\bm{\phi}_{j}(i)=(\phi_{ij,1},\cdots,\phi_{ij,j})^{*}\in\mathbb{R}^{j}.$ It is easy to see from (1) of Assumption \arabicsection.\arabicthm and Yule-Walker’s equation that

\bm{\phi}_{j}(i)=\Gamma_{j}^{-1}(i)\bm{\nu}_{j}(i),

(A.\arabicequation)

where we denoted that

\Gamma_{j}(i)=\operatorname{Cov}\left(\bm{x}_{j}(i),\bm{x}_{j}(i)\right)\in% \mathbb{R}^{j\times j},\ \bm{\nu}_{j}(i)=\operatorname{Cov}\left(\bm{x}_{j}(i)% ,x_{i}\right)\in\mathbb{R}^{j},

where $\bm{x}_{j}(i)=(x_{i-1},\cdots,x_{i-j})\in\mathbb{R}^{j}.$ For the rest of the proof, we can follow lines of that of Theorem 2.4 of Ding & Zhou (2023) verbatim. More specifically, according to (A.\arabicequation), we find that in order to study $\rho_{i,j}\equiv\phi_{ij,j},$ it suffices to control the entries of the $j$ th row of $\Gamma_{j}^{-1}(i)$ and all the entries of $\bm{\nu}_{j}(i).$ For $\Gamma_{j}^{-1}(i),$ we denote the $j\times j$ symmetric banded matrix $\Gamma_{j}^{s}(i)$ as

(\Gamma_{j}^{s}(i))_{kl}=\begin{cases}(\Gamma_{j}(i))_{kl},&|k-l|\leq\frac{j}{% K\log j};\\ 0,&\text{otherwise}.\end{cases}

Here $K>0$ is some large constant. Using (1) and (2) of Assumption \arabicsection.\arabicthm, by a discussion similar to equation (D.7) of Ding & Zhou (2023), we find that for some constant $C>0$

|\rho_{i,j}-\phi_{ij,j}^{s}|\leq Cj^{1-\tau}(K\log j)^{\tau-1},

where $\bm{\phi}_{j}(i)^{s}=(\phi_{ij,1}^{s},\cdots,\phi_{ij,j}^{s})^{*}\in\mathbb{R}% ^{j}$ is defined via $\bm{\phi}_{j}^{s}(i)=(\Gamma_{j}(i)^{s})^{-1}\bm{\nu}_{j}(i).$ Moreover, using (1) of Assumption \arabicsection.\arabicthm, according to the discussions below equation (D.8) of Ding & Zhou (2023), we find that for some constant $C_{1}>0,$

|\phi_{ij,j}^{s}|\leq C_{1}(j/(\log j+1))^{-\tau+1}.

Combining the above controls, we complete the proof of (\arabicsection.\arabicequation).

For the second part of the proof, for the smoothness that $\rho_{j}(t)\in C^{d}([0,1])$ follows from Lemma 3.1 of Ding & Zhou (2020). To see (\arabicsection.\arabicequation), by (A.\arabicequation) and (\arabicsection.\arabicequation), using Cauchy-Schwarz inequality, we have that

\displaystyle\left|\bm{\phi}_{j}(i)-\bm{\phi}_{j}(i/n)\right|\leq\left\|\Gamma% _{j}^{-1}(i)\right\|\left|\bm{\nu}_{j}(i)-\bm{\nu}_{j}(i/n)\right|+\left\|% \Gamma_{j}^{-1}(i)-\Gamma_{j}^{-1}(i/n)\right\|\left|\bm{\nu}_{j}(i/n)\right|.

(A.\arabicequation)

The rest of the proof follows from the exact reasoning as the arguments of the proof of Theorem 2.11 of Ding & Zhou (2023). In particular, using (1)-(3) of Assumption \arabicsection.\arabicthm, by an argument similar to the equations between (D.31) and (D.32) of Ding & Zhou (2023), we find that the first term of the right-hand side of (A.\arabicequation) can be bounded by $\mathrm{O}(j^{1.5}/n)$ and the second term can be bounded by $\mathrm{O}(j^{2}/n).$ This completes our proof.

A.\arabicsubsection Proofs of Section \arabicsection

In this subsection, we prove the main results in Section \arabicsection.

{Proof}

[Proof of Theorem \arabicsection.\arabicthm] The proof is similar to that of Theorem 3.2 of Ding & Zhou (2021). Due to similarity, we only sketch the key points. Similar to (\arabicsection.\arabicequation), we set $\widehat{\bm{\phi}}_{j}(t)$ as the sieve estimator for $\bm{\phi}_{j}(t);$ that is

\widehat{\bm{\phi}}_{j}(t):=\widehat{\bm{\beta}}_{j}^{*}\mathbb{B}_{j}(t),\ % \mathbb{B}_{j}(t)=\sum_{l=1}^{j}\mathbb{B}_{j,l}(t).

(A.\arabicequation)

Recall (\arabicsection.\arabicequation). For $1\leq l\leq j,$ denote

\bm{\phi}_{j,l}^{(c)}(t)=\sum_{k=1}^{c}a_{jk,l}\alpha_{k}(t),

and $\bm{\phi}_{j}^{(c)}=\left(\bm{\phi}_{j,1}^{(c)},\cdots,\bm{\phi}_{j,j}^{(c)}% \right)^{*}\in\mathbb{R}^{j}.$ Recall (\arabicsection.\arabicequation). The starting point is the following decomposition

\|\widehat{\bm{\phi}}_{j}(t)-\bm{\phi}_{j,j}(t)\|\leq\|\widehat{\bm{\phi}}_{j}% (t)-\bm{\phi}^{(c)}_{j,j}(t)\|+|\bm{\phi}_{j,j}(t)-\bm{\phi}_{j,j}^{(c)}(t)|.

Since the second term of the right-hand side of the above equation can be bounded by $\mathrm{O}\left(c^{-d}\right)$ using (\arabicsection.\arabicequation), it suffices to control the first term. According to (A.\arabicequation), by Cauchy-Schwarz inequality, we see that

\|\widehat{\bm{\phi}}_{j,j}(t)-\bm{\phi}^{(c)}_{j,j}(t)\|\leq\sum_{l=1}^{j}\|% \widehat{\bm{\phi}}_{j,l}(t)-\bm{\phi}^{(c)}_{j,l}(t)\|\leq\sqrt{j}\zeta_{c}\|% \bm{\beta}_{j}-\widehat{\bm{\beta}}_{j}\|.

(A.\arabicequation)

Moreover, for the OLS estimation, we have that

\bm{\beta}_{j}-\widehat{\bm{\beta}}_{j}=n(Y_{j}^{*}Y_{j})^{-1}\frac{Y_{j}^{*}% \bm{\epsilon}_{j}}{n},

(A.\arabicequation)

where $\bm{\epsilon}_{j}=(\epsilon_{j+1,j}+\mathfrak{r},\cdots,\epsilon_{n,j}+% \mathfrak{r})^{*}\in\mathbb{R}^{n-j}$ with $\mathfrak{r}$ representing the error term in (\arabicsection.\arabicequation). The rest of the discussions follow lines of that of Theorem 3.2 of Ding & Zhou (2021).

To control the right-hand side of (A.\arabicequation), first, by an argument similar to (A.9) of Ding & Zhou (2021), for $\Sigma^{(j)}$ defined in (\arabicsection.\arabicequation), we find that

\left\|\frac{1}{n}Y_{j}^{*}Y_{j}-\Sigma^{(j)}\right\|=\mathrm{O}_{\mathbb{P}}% \left(jc\left(\frac{\xi_{c}^{2}}{\sqrt{n}}+\frac{\xi_{c}^{2}n^{\frac{2}{\tau+1% }}}{n}\right)\right).

Combining the above discussion, together with Assumption \arabicsection.\arabicthm and the assumption of (\arabicsection.\arabicequation), we see that

\left\|n(Y_{j}^{*}Y_{j})^{-1}\right\|=\mathrm{O}_{\mathbb{P}}(1).

(A.\arabicequation)

Second, according to a discussion similar to (A.12) of Ding & Zhou (2021) and the fact that $\mathfrak{r}=\mathrm{O}_{\mathbb{P}}(j^{2}/n+jc^{-d})$ , we have that

\left\|\frac{Y_{j}^{*}\bm{\epsilon}_{j}}{n}\right\|=\mathrm{O}_{\mathbb{P}}% \left(\xi_{c}\left(1+\frac{j^{2}}{n}+jc^{-d}\right)\sqrt{\frac{jc}{n}}\right).

(A.\arabicequation)

We point out that the error rate $j^{2}/n$ is faster than the $j^{2.5}/n$ in (A.12) of Ding & Zhou (2021) or (2.23) and (2.24) of Ding & Zhou (2023) because of the assumption that $\{x_{i}\}$ is a mean-zero time series. Inserting the above two bounds into (A.\arabicequation), we conclude that

\|\bm{\beta}_{j}-\widehat{\bm{\beta}}_{j}\|=\mathrm{O}_{\mathbb{P}}\left(\xi_{% c}\left(1+\frac{j^{2}}{n}+jc^{-d}\right)\sqrt{\frac{jc}{n}}\right).

Combining (A.\arabicequation), we can conclude our proof.

{Proof}

[Proof of Theorem \arabicsection.\arabicthm] The proof follows from lines of those of Proposition 3.7 of Ding & Zhou (2023) verbatim. We omit the details.

{Proof}

[Proof of Theorem \arabicsection.\arabicthm] The proof of part 1 follows from lines of those of Proposition 3.7 of Ding & Zhou (2023) verbatim. We omit the details.

For part 2, we start with (\arabicsection.\arabicequation). Before proceeding to the actual proof, we first explore the relation between $\rho_{k}(t)$ and $\phi_{\mathsf{h},k}(t)$ for $1\leq k\leq\mathsf{h}$ under the null hypothesis of (\arabicsection.\arabicequation). First of all, under (\arabicsection.\arabicequation), according to Theorem \arabicsection.\arabicthm, we can conclude that

\rho_{i,\mathsf{h}}\equiv\phi_{i\mathsf{h},\mathsf{h}}=\mathrm{O}\left(\min% \left\{\frac{\mathsf{h}^{2}}{n},\left(\frac{\log\mathsf{h}+1}{\mathsf{h}}% \right)^{\tau-1}\right\}\right).

(A.\arabicequation)

Combining with an argument similar to Theorem 2.11 of Ding & Zhou (2023) and the smoothness of $\phi_{\mathsf{h},\mathsf{h}}(t)$ , we can conclude that for all $0\leq t\leq 1$

\phi_{\mathsf{h},\mathsf{h}}(t)=\mathrm{O}\left(\frac{\mathsf{h}^{2}}{n}\right).

(A.\arabicequation)

Then with an argument similar to Theorem 2.11 of Ding & Zhou (2023), in light of the representation (\arabicsection.\arabicequation), we see that

x_{i}=\sum_{j=1}^{\mathsf{h}-1}\phi_{\mathsf{h},j}\left(\frac{i}{n}\right)x_{i% -j}+\epsilon_{i,\mathsf{h}}+\mathrm{O}_{\mathbb{P}}\left(\frac{\mathsf{h}^{2}}% {n}\right).

Therefore, we can repeat (A.\arabicequation) and (A.\arabicequation) to get that

\phi_{\mathsf{h},k}(t)=\mathrm{O}\left(\frac{k^{2}}{n}\right),\ 1\leq k\leq% \mathsf{h},

(A.\arabicequation)

and conclude that

x_{i}=\epsilon_{i,\mathsf{h}}+\mathrm{O}_{\mathbb{P}}\left(\frac{\mathsf{h}^{3% }}{n}\right).

Using our assumptions on $\mathsf{h}$ , we find that the error term of the above approximation is always negligible with high probability. With the above preparation, we proceed to the proof. By definition, we have that

T_{\text{BP}}-T_{2}=\sum_{k=1}^{\mathsf{h}}\int_{0}^{1}\left(\widehat{\rho}_{k% }^{2}(t)-\widehat{\phi}_{\mathsf{h},k}^{2}(t)\right)\mathrm{d}t.

By a discussion similar to Theorem \arabicsection.\arabicthm, we have that (see Theorem 3.3 of Ding & Zhou (2020)) for $1\leq k\leq\mathsf{h}$

\phi_{\mathsf{h},k}(t)=\widehat{\phi}_{\mathsf{h},k}(t)+\mathrm{O}_{\mathbb{P}% }\left(\Psi(\mathsf{h},c)\right),

where we recall (\arabicsection.\arabicequation). Combining the above controls with Theorem \arabicsection.\arabicthm, under the null hypothesis (\arabicsection.\arabicequation), we see that

\widehat{\rho}_{k}^{2}(t)=\mathrm{O}_{\mathbb{P}}\left(\Psi^{2}(k,c)\right),\ % \widehat{\phi}_{\mathsf{h},k}^{2}(t)=\mathrm{O}_{\mathbb{P}}\left(\frac{k^{4}}% {n^{2}}+\Psi^{2}(\mathsf{h},c)\right).

Consequently, we have that

T_{\text{BP}}-T_{2}=\mathrm{O}_{\mathbb{P}}\left(\mathsf{h}\Psi^{2}(\mathsf{h}% ,c)+\frac{\mathsf{h}^{5}}{n^{2}}+\sum_{k=1}^{\mathsf{h}}\Psi^{2}(k,c)\right).

Then, for (\arabicsection.\arabicequation), using the definition of $g_{k}$ in (\arabicsection.\arabicequation) and the fact that $g_{2}\asymp\sqrt{\mathsf{h}c},$ the proof follows from (\arabicsection.\arabicequation), (\arabicsection.\arabicequation) and the assumption of (\arabicsection.\arabicequation).

Finally, for part 3, by a discussion similar to (A.\arabicequation) and (A.\arabicequation), we find that (\arabicsection.\arabicequation) yields that

\sum_{k=1}^{\mathsf{h}}\int_{0}^{1}\phi_{\mathsf{h},k}^{2}(t)\mathrm{d}t>C_{% \alpha}\frac{\sqrt{\mathsf{h}c}}{n}.

Then the proof follows from lines of those of Proposition 3.7 of Ding & Zhou (2023) verbatim.

{Proof}

[Proof of Corollary \arabicsection.\arabicthm] The proof follows directly from the proof of Theorem 3.10 of Ding & Zhou (2023) by setting $b_{*}=\mathsf{h}$ therein.

A.\arabicsubsection Additional assumptions

In this section, we collect some more assumptions and provide more discussions on these technical assumptions. The following assumption will be needed in our proof. It is a mild assumption and can be easily satisfied by many time series. We refer the readers to Section C.1 of Ding & Zhou (2023) for more details.

Assumption A.\arabicthm.

We assume the following assumptions hold true

\arabicenumi.

Suppose $\tau>4.5$ in Assumption \arabicsection.\arabicthm. Moreover, we assume that $c$ is of the form (\arabicsection.\arabicequation) and satisfies that for large constant $C>0$

\frac{C}{\tau}+a<1,\ \operatorname{and}\ da>2.

\arabicenumi.

We assume that the derivatives of $\gamma(t,j)$ decay with $j$ as follows

\sup_{t\in[0,1]}\sum_{j=0}^{\infty}|\gamma^{(d)}(t,j)|<\infty,

where $\gamma^{(d)}(t,j)$ is the d $th$ derivative of $\gamma(t,j)$ with respect to $t$ .

\arabicenumi.

There exist constants $\omega_{1},\omega_{2}\geq 0,$ for some constant $C>0,$ we have

\sup_{t}|\nabla\mathbf{B}(t)|\leq Cn^{\omega_{1}}c^{\omega_{2}}.

B Additional remarks and examples

In this section, we provide several remarks and examples. The following remark provides more discussions on Theorem \arabicsection.\arabicthm.

Remark B.\arabicthm.

Theorem \arabicsection.\arabicthm is established for locally stationary time series as in Definition \arabicsection.\arabicthm where an exact cut-off is not available in general. An exception is the locally stationary AR(p) process as in Zhou (2013b), where

x_{i}=\phi_{0}(i/n)+\sum_{j=1}^{p}\phi_{j}(i/n)x_{i-j}+\epsilon_{i},

(B.\arabicequation)

where $\{\epsilon_{i}\}$ is a time-varying white noise process. In this setting, it is clear that

\rho_{i,j}=\begin{cases}\phi_{j}(i/n),&j=p;\\ 0,&j>p,\end{cases}

so that $(\ref{eq_decayrate})$ and (\arabicsection.\arabicequation) holds trivially once $p$ is fixed or divergent slowly. In fact, as proved in Theorem 2.11 of Ding & Zhou (2023), any locally stationary time series satisfying Definition \arabicsection.\arabicthm and Assumption \arabicsection.\arabicthm can be always well approximated by a time series in the form of (B.\arabicequation) with slowly diverging $p.$ In this regard, our results can be used to provide an order for the AR approximation; see Remark B.\arabicthm below.

The following remark provides more explanations on the conditions in Assumption \arabicsection.\arabicthm.

Remark B.\arabicthm.

The conditions (1)–(3) in Assumption \arabicsection.\arabicthm are mild and commonly used in the literature. First, (1) is introduced to avoid the erratic behavior of the time series and frequently used in the statistics literature involving the covariance and precision matrix estimation Cai et al. (2016); Chen et al. (2013); Ding & Zhou (2020, 2023); Yuan (2010). Moreover, as proved in (Ding & Zhou, 2023, Proposition 2.9), it is equivalent to the uniform positiveness of the local spectral density function of $\{x_{i}\}.$ Second, (2) imposes the condition that the temporal structure of $\{x_{i}\}$ decays polynomially fast. This amounts to a short range dependent requirement for $\{x_{i}\}$ when $\tau>1.$ Analogous results can be easily obtained for the exponentially decay setting where

\max_{k,n}\left|\operatorname{Cov}(x_{k,n},x_{k+r,n})\right|\leq Ca^{|r|},\ % \text{for some}\ 0<a<1.

Third, (3) requires that the autocovariance functions of $\{x_{i}\}$ are smooth so that its PACFs can be estimated consistently. It is commonly used in the literature of locally stationary time series Dahlhaus (2012); Dahlhaus et al. (2019); Ding & Zhou (2020, 2023); Zhou & Wu (2009).

The remark below provides some insights on how to use our results to estimate the order of a locally stationary AR process.

Remark B.\arabicthm.

We discuss how to generalize the use of PACFs of stationary AR process to locally stationary AR process. For definiteness, we focus on the following time-varying AR( $p$ ) process which has been used in Zhou (2013b); Dahlhaus et al. (2019)

x_{i}=\sum_{j=1}^{p}\phi_{j}(i/n)x_{i-j}+\epsilon_{i},

(B.\arabicequation)

where $\{\epsilon_{i}\}$ is some locally-stationary white noise process and $\phi_{j}(\cdot),j=0,1,2,\cdots,p,$ are some smooth functions on $[0,1]$ . Before estimating these time-varying coefficients, we first need to provide an estimator for $p.$ Here $p$ is allowed to diverge with $n.$ Based on the results of Theorem \arabicsection.\arabicthm, inspired by the ideas in Ding & Yang (2022); Ding et al. (2023), we can propose a sequential test estimate based on the following hypothesis testing problem

\mathbf{H}_{0}:p=p_{0}\ \text{vs}\ \mathbf{H}_{a}:p_{0}<p\leq p_{*},

(B.\arabicequation)

where $p_{0}$ is some pre-given integer representing our belief of the true value of $p$ and $p_{*}$ is some large integer that is interpreted as the maximum possible order the model can have. In light of (\arabicsection.\arabicequation), we can use the following estimate

\widetilde{p}=\max\{1\leq j\leq p_{*}:\ \mathbf{H}_{a1}\ \text{in}\ (\ref{eq_% firstnull})\ \text{is accepted}\}.

(B.\arabicequation)

The following remark is related to the hypothesis testing (\arabicsection.\arabicequation).

Remark B.\arabicthm.

Two remarks are in order. First, in classic stationary time series analysis, in addition to Box-Pierce test, people also use the Ljung–Box (LB) test Ljung & Box (1978). Moreover, the BP and LB tests are asymptotic equivalent and follow the Chi-squared distribution with the same degree of freedom $\mathsf{h}$ . In this regard, we can also modify the LB test using

\mathcal{Q}_{\text{MLB}}:=n(n+2)\sum_{k=1}^{\mathsf{h}}\frac{\int_{0}^{1}% \widehat{\rho}_{k}(t)^{2}}{n-k}.

Moreover, we can study such a modified statistic as Theorem \arabicsection.\arabicthm. Due to the asymptotic equivalence, we omit further details. Second, (\arabicsection.\arabicequation) is frequently used for model diagnostics. In this regard, it provides an alternative approach to choose the order of AR approximations by checking whether the residuals follow white noise after fitting some AR models.

Finally, we provide two frequently-used models of locally stationary time series in the literature and explain how Definition \arabicsection.\arabicthm and Assumption \arabicsection.\arabicthm can be easily satisfied.

Example B.\arabicthm.

We shall first consider the locally stationary time series model in Zhou & Wu (2009, 2010) using a physical representation so that

x_{i,n}=G_{n}(\frac{i}{n},\mathcal{F}_{i}),\ i=1,2,\cdots,n,

(B.\arabicequation)

where $\mathcal{F}_{i}=(\cdots,\eta_{i-1},\eta_{i})$ and $\eta_{i},\ i\in\mathbb{Z}$ are i.i.d centered random variables, and $G_{n}:[0,1]\times\mathbb{R}^{\infty}\rightarrow\mathbb{R}$ is a measurable function such that $\xi_{i,n}(t):=G_{n}(t,\mathcal{F}_{i})$ is a properly defined random variable for all $t\in[0,1].$ In (B.\arabicequation), by allowing the data generating mechanism $G_{n}$ depending on the time index $t$ in such a way that $G_{n}(t,\mathcal{F}_{i})$ changes smoothly with respect to $t$ , one has local stationarity in the sense that the subsequence $\{x_{i,n},...,x_{i+j-1,n}\}$ is approximately stationary if its length $j$ is sufficiently small compared to $n$ . Moreover, they quantify the temporal decay using the physical dependence measure for (B.\arabicequation) as follows

\delta(j,q):=\sup_{t\in[0,1]}||G_{n}(t,\mathcal{F}_{0})-G_{n}(t,\mathcal{F}_{0% ,j})||_{q}.

(B.\arabicequation)

Moreover, the following assumptions are needed to ensure local stationarity.

Assumption B.\arabicthm.

$G_{n}(\cdot,\cdot)$ defined in (B.\arabicequation) satisfies the property of stochastic Lipschitz continuity, i.e., for some $q>2$ and $C>0,$

\left|\left|G_{n}(t_{1},\mathcal{F}_{i})-G_{n}(t_{2},\mathcal{F}_{i})\right|% \right|_{q}\leq C|t_{1}-t_{2}|,

(B.\arabicequation)

where $t_{1},t_{2}\in[0,1].$ Furthermore,

\sup_{t\in[0,1]}\max_{1\leq i\leq n}||G_{n}(t,\mathcal{F}_{i})||_{q}<\infty.

(B.\arabicequation)

It can be shown that time series $\{x_{i,n}\}$ with physical representation (B.\arabicequation) and Assumption B.\arabicthm satisfies Definition \arabicsection.\arabicthm. In particular, for each fixed $t\in[0,1],$ $\gamma(t,j)$ in Definition \arabicsection.\arabicthm can be found easily using the following

\gamma(t,j)=\operatorname{Cov}(G_{n}(t,\mathcal{F}_{0}),G_{n}(t,\mathcal{F}_{j% })).

(B.\arabicequation)

Note that the assumptions (B.\arabicequation) and (B.\arabicequation) ensure that $\gamma(t,j)$ is Lipschiz continuous in $t$ . Moreover, for each fixed $t,$ $\gamma(t,\cdot)$ is the autocovariance function of $\{G_{n}(t,\cdot)\},$ which is a stationary process.

The physical representation form (B.\arabicequation) includes many commonly used locally stationary time series models. For example, let $\{z_{i}\}$ be zero-mean i.i.d. random variables (or a white noise) with variance $\sigma^{2}$ . We also assume $a_{j,n}(\cdot),j=0,1,\cdots$ be $C^{d}([0,1])$ functions such that

G_{n}(t,\mathcal{F}_{i})=\sum_{k=0}^{\infty}a_{k,n}(t)z_{i-k}.

(B.\arabicequation)

(B.\arabicequation) is a locally stationary linear process. It is easy to see that (2) and (3) of Assumption \arabicsection.\arabicthm will be satisfied if $\sup_{t\in[0,1]}|a_{j,n}(t)|\leq Cj^{-\tau},\ j\geq 1;\ \sum_{j=0}^{\infty}% \sup_{t\in[0,1]}|a_{j,n}^{\prime}(t)|<\infty,$ and

\sup_{t\in[0,1]}|a_{j,n}^{(d)}(t)|\leq Cj^{-\tau},\ j\geq 1.

(B.\arabicequation)

Furthermore, we note that the local spectral density function of (B.\arabicequation) can be written as $f(t,w)=\sigma^{2}|\psi(t,e^{-\mathrm{i}j\omega})|^{2},$ where $\psi(\cdot,\cdot)$ is defined such that $G_{n}(t,\mathcal{F}_{i})=\psi(t,B)z_{i}$ with $B$ being the backshift operator. As discussed in Remark B.\arabicthm, (1) of Assumption B.\arabicequation will be satisfied if $|\psi(t,e^{-\mathrm{i}j\omega})|^{2}\geq\kappa$ for all $t$ and $\omega,$ where $\kappa>0$ is some universal constant. For more examples of locally stationary time series in the form of (B.\arabicequation) especially nonlinear time series, we refer the readers to Wu (2005), (Ding & Zhou, 2020, Section 2.1) , (Dahlhaus et al., 2019, Example 2.2 and Proposition 4.4), (Karmakar et al., 2022, Proposition E.6) and Ding & Zhou (2021); Karmakar et al. (2022); Mayer et al. (2020). Especially, the time-varying AR and ARCH models can be written into (B.\arabicequation) asymptotically Ding & Zhou (2023), and Assumptions \arabicsection.\arabicthm and B.\arabicthm can be easily satisfied under mild assumptions. We refer the readers to the aforementioned references for more details.

For a second example, note that in Dahlhaus & Rao (2006); Dette et al. (2011); Vogt (2012), the locally stationary time series is defined as follows (see Definition 2.1 of Vogt (2012)). $\{x_{i,n}\}$ is locally stationary time series if for each scaled time point $u\in[0,1],$ there exists a strictly stationary process $\{h_{i,n}(u)\}$ such that

|x_{i,n}-h_{i,n}(u)|\leq\left(|t_{i}-u|+\frac{1}{n}\right)U_{i,n}(u),\ \text{a% .s},

(B.\arabicequation)

where $U_{i,n}(u)\in L^{q}([0,1])$ for some $q>0.$ By similar arguments as those of model (B.\arabicequation) Ding & Zhou (2023), Definition \arabicsection.\arabicthm as well as assumptions of this subsection can be verified for (B.\arabicequation), especially (B.\arabicequation) implies (\arabicsection.\arabicequation).

C Tuning parameters selection

In this section, we explain how to choose the tuning parameters associated with our proposed methodology.

First, we discuss how to choose the tuning parameters $c$ and $m$ used in Algorithm 1. We use a data-driven procedure proposed in Bishop (2013) to choose $c.$ For a given integer $l,$ say $l=\lfloor 3\log_{2}n\rfloor,$ we divide the time series into two parts: the training part $\{x_{i}\}_{i=1}^{n-l}$ and the validation part $\{x_{i}\}_{i=n-l+1}^{n}.$ With some preliminary initial value $c$ , we propose a sequence of candidates $c_{j},\ j=1,2,\cdots,v,$ in an appropriate neighborhood of $c$ where $v$ is some given integer. For each of the choices $c_{j},$ we fit a time-varying AR( $h_{2}$ ) model as in (\arabicsection.\arabicequation) with $c_{j}$ sieve basis expansion using the training data set. Then using the fitted model, we forecast the time series in the validation part of the time series. Let $\widehat{x}_{n-l+1,j},\cdots,\widehat{x}_{n,j}$ be the forecast of $x_{n-l+1},...,x_{n},$ respectively using the parameter $c_{j}$ . Then we choose the parameter $c_{j_{0}}$ with the minimum sample MSE of forecast, i.e.,

{j_{0}}:=\operatorname*{argmin}_{(j:1\leq j\leq v)}\frac{1}{l}\sum_{k=n-l+1}^{% n}(x_{k}-\widehat{x}_{k,j})^{2}.

To choose an $m$ for practical implementation, in Zhou (2013a), the author used the minimum volatility (MV) method to choose the window size $m$ for the scalar covariance function. The MV method does not depend on the specific form of the underlying time series dependence structure and hence is robust to misspecification of the latter structure Politis et al. (1999). The MV method utilizes the fact that the covariance structure of $\widehat{\Pi}$ becomes stable when the block size $m$ is in an appropriate range, where $\widehat{\Pi}=E[\widehat{\Phi}\widehat{\Phi}^{*}|(x_{1},\cdots,x_{n})]$ is defined as

\widehat{\Pi}:=\frac{1}{(n-m-h_{2}+1)m}\sum_{i=h_{2}+1}^{n-m}\Big{[}\Big{(}% \sum_{j=i}^{i+m}\widehat{\bm{w}}_{h_{2},j}\Big{)}\otimes\Big{(}\mathbf{B}(% \frac{i}{n})\Big{)}\Big{]}\times\Big{[}\Big{(}\sum_{j=i}^{i+m}\widehat{\bm{w}}% _{h_{2},j}\Big{)}\otimes\Big{(}\mathbf{B}(\frac{i}{n})\Big{)}\Big{]}^{*}.

(C.\arabicequation)

Therefore, it desires to minimize the standard errors of the latter covariance structure in a suitable range of candidate $m$ ’s. In detail, for a give large value $m_{n_{0}}$ and a neighborhood control parameter $h_{0}>0,$ we can choose a sequence of window sizes $m_{-h_{0}+1}<\cdots<m_{1}<m_{2}<\cdots<m_{n_{0}}<\cdots<m_{n_{0}+h_{0}}$ and obtain $\widehat{\Pi}_{m_{j}}$ by replacing $m$ with $m_{j}$ in (C.\arabicequation), $j=-h_{0}+1,2,\cdots,n_{0}+h_{0}.$ For each $m_{j},j=1,2,\cdots,m_{n_{0}},$ we calculate the matrix norm error of $\widehat{\Omega}_{m_{j}}$ in the $h_{0}$ -neighborhood, i.e.,

\mathsf{se}(m_{j}):=\mathsf{se}(\{\widehat{\Pi}_{m_{j+k}}\}_{k=-h_{0}}^{h_{0}}% )=\left[\frac{1}{2h_{0}}\sum_{k=-h_{0}}^{h_{0}}\|\overline{\widehat{\Pi}}_{m_{% j}}-\widehat{\Pi}_{m_{j}+k}\|^{2}\right]^{1/2},

where $\overline{\widehat{\Pi}}_{m_{j}}=\sum_{k=-h_{0}}^{h_{0}}\widehat{\Pi}_{m_{j}+k% }/(2h_{0}+1).$ Therefore, we choose the estimate of $m$ using

\widehat{m}:=\operatorname*{argmin}_{m_{1}\leq m\leq m_{n_{0}}}\mathsf{se}(m).

Note that in Zhou (2013a) the author used $h_{0}=3$ and we also adopt this choice in the current paper.

Second, we discuss how to choose a large value of $\mathsf{h}$ to construct the statistic $T_{2}$ in (\arabicsection.\arabicequation). Theoretically, the lower bound for the order of $\mathsf{h}$ is given by $j^{*}$ as in (\arabicsection.\arabicequation), while the upper bound is provided in the assumption of (\arabicsection.\arabicequation). The lower bound yields that for all $j\geq\mathsf{h},$ $\sup_{t}\rho_{j}(t)=\mathrm{o}(n^{-1/2}),$ while the upper bound guarantees that the error term is negligible. To balance these two conditions, for practical implementation, we use the following value

\mathsf{h}:=\min\left\{1\leq j\leq\mathsf{h}^{*}:\ \mathbf{H}_{01}\ \text{in}% \ (\ref{eq_firstnull})\ \text{is accepted}\right\},

where $\mathsf{h}^{*}>0$ is some pre-given large value (say, $\mathsf{h}^{*}=50$ ). We emphasize that this can be easily done using our $\mathtt{R}$ package $\mathtt{Sie2nts}$ by generating a plot as in Figure \arabicfigure.

D Additional simulation results

In this section, we provide additional numerical simulation results.

D.\arabicsubsection More results on other types of models

In this section, we conduct more numerical simulations using both stationary and non-stationary MA(1) models. For some constant $\delta\in[0,0.5],$ we consider the stationary MA(1) process

x_{i}=\epsilon_{i}+\delta\epsilon_{i-1},

(D.\arabicequation)

and the locally stationary MA(1) process

x_{i}=\epsilon_{i}+\delta\sin(2\pi i/n)\epsilon_{i-1},

(D.\arabicequation)

where $\epsilon_{i},1\leq i\leq n,$ are i.i.d. standard Gaussian random variables. Note that when $\delta=0$ in (D.\arabicequation) and (D.\arabicequation), they both reduce to the standard white noise. Since (D.\arabicequation) and (D.\arabicequation) are essentially AR( $\infty$ ) models, one can follow the discussions of Section 3.3 of Shumway & Stoffer (2017) to calculate the true PACFs.

For the estimation of the PACFs, in Figure A, we provide the plots of the PACFs of the first 10 lags for both (D.\arabicequation) and (D.\arabicequation). We can see that our estimates are reasonably accurate. Regarding the inference of the PACFs, for the purpose of definiteness, we focus on the white noise test (\arabicsection.\arabicequation) where $\mathbf{H}_{0}$ corresponds to $\delta=0$ in (D.\arabicequation) and (D.\arabicequation) and $\mathbf{H}_{a}$ corresponds to an MA(1) alternative that $\delta>0.$ Under $\mathbf{H}_{0},$ (D.\arabicequation) and (D.\arabicequation) are essentially the same model. For $n=600,$ under the type I error rate $\alpha=0.05,$ the simulated type I error rates are $0.047,0.051,0.056$ for the Fourier, Legendre and Daubechies-9 basis functions, respectively based on 1,000 repetitions. This shows the accuracy of our test. To examine the power, in Figure B, we report how the simulated power changes when $\delta$ deviates away from zero. We can conclude that our proposed test is reasonably powerful once the alternative deviates from the null.

D.\arabicsubsection Comparison with Killick et al. (2020) on estimating the PACFs

In this section, we compare our method with the ones proposed in Killick et al. (2020) in terms of the estimation of the PACFs using the mean integrated squared error (MISE). To implement Killick et al. (2020), we use the $\mathtt{R}$ package $\mathtt{lpacf}$ which is developed by the authors of Killick et al. (2020).

For definiteness, we consider the AR type models (\arabicsection.\arabicequation) and (\arabicsection.\arabicequation) with $\delta_{1}=0.5$ and $\delta_{2}=0.$ Consequently, for model (\arabicsection.\arabicequation), the true PACFs are

\rho_{j}(t)\equiv\begin{cases}0.5&j=1\\ 0&j\geq 2\end{cases},

and for model (\arabicsection.\arabicequation), the true PACFs are

\rho_{j}(t)=\begin{cases}0.5\sin(2\pi t)&j=1\\ 0&j\geq 2\end{cases}.

Let $\widehat{\rho}_{j}(t)$ be some estimator for $\rho_{j}(t),$ the MISE is defined as

\mathsf{MISE}(j)=\int_{0}^{1}(\widehat{\rho}_{j}(t)-\rho_{j}(t))^{2}\mathrm{d}t.

In the actual calculations, MISE can be well approximated by an Riemann summation. In Table A, we record the $\mathsf{MISE}(j),j=1,2,3,4,$ for our proposed method (denoted as Proposed) and the two methods in Killick et al. (2020) (the wavelet-based method is denoted as Lpacf-I and the Epanechnikov windowed method is denoted as Lpacf-II). We can conclude that our proposed method has better finite sample performance than Killick et al. (2020) which are known to have worse performance near the boundaries.

Methods/Lags	$j=1$	$j=2$	$j=3$	$j=4$
	Model (\arabicsection.\arabicequation)
Proposed	$8\times 10^{-4}$	$7\times 10^{-4}$	$3\times 10^{-3}$	$1\times 10^{-3}$
Lpacf-I	$9\times 10^{-3}$	$7\times 10^{-3}$	$9\times 10^{-3}$	$7\times 10^{-3}$
Lpacf-II	0.042	0.054	0.043	0.051
	Model (\arabicsection.\arabicequation)
Proposed	$9.8\times 10^{-3}$	$1.9\times 10^{-3}$	$2.5\times 10^{-3}$	$2.7\times 10^{-3}$
Lpacf-I	0.016	$7.5\times 10^{-3}$	$3.5\times 10^{-3}$	$4\times 10^{-3}$
Lpacf-II	0.049	0.051	0.042	0.055

Table A: Comparison of the accuracy of the estimation of PACFs using MISE. Here Proposed stands for our estimator (\arabicsection.\arabicequation) and Lpacf-I and Lpacf-II are the two methods introduced in Killick et al. (2020). Our proposed method can be implemented using our package

\mathtt{Sie2nts}

and Lpacf-I/II can be implemented using the package

\mathtt{lpacf}.

Here

n=1024.

References

Bishop (2013) Bishop, C. (2013). Pattern Recognition and Machine Learning. Information science and statistics. Springer.
Box & Pierce (1970) Box, G. E. & Pierce, D. A. (1970). Distribution of residual autocorrelations in autoregressive-integrated moving average time series models. Journal of the American statistical Association 65, 1509–1526.
Brockwell & Davis (1987) Brockwell, P. & Davis, R. (1987). Time series: Theory and Methods. Springer-Verlag.
Cai et al. (2016) Cai, T. T., Liu, W. & Zhou, H. H. (2016). Estimating sparse precision matrix: Optimal rates of convergence and adaptive estimation. Ann. Statist. 44, 455–488.
Chen (2007) Chen, X. (2007). Large sample sieve estimation of semi-nonparametric models. Handbook of econometrics 6, 5549–5632.
Chen et al. (2013) Chen, X., Xu, M. & Wu, W. B. (2013). Covariance and precision matrix estimation for high-dimensional time series. Ann. Statist. 41, 2994–3021.
Dahlhaus (2012) Dahlhaus, R. (2012). Locally stationary processes. In Handbook of statistics, vol. 30. Elsevier, pp. 351–413.
Dahlhaus & Rao (2006) Dahlhaus, R. & Rao, S. S. (2006). Statistical inference for time-varying arch processes. The Annals of Statistics 34, 1075–1114.
Dahlhaus et al. (2019) Dahlhaus, R., Richter, S. & Wu, W. B. (2019). Towards a general theory for nonlinear locally stationary processes. Bernoulli 25, 1013–1044.
Dette et al. (2011) Dette, H., Preubb, P. & Vetter, M. (2011). A measure of stationarity in locally stationary processes with applications to testing. J. Am. Stat. Assoc. 106, 1113–1124.
Dette & Wu (2020) Dette, H. & Wu, W. (2020). Prediction in locally stationary time series. Journal of Business & Economic Statistics , 1–12.
Ding et al. (2023) Ding, X., Xie, J., Yu, L. & Zhou, W. (2023). Extreme eigenvalues of sample covariance matrices under generalized elliptical models with applications. arXiv preprint arXiv:2303.03532 .
Ding & Yang (2022) Ding, X. & Yang, F. (2022). Tracy-Widom distribution for heterogeneous gram matrices with applications in signal detection. IEEE Transactions on Information Theory 68, 6682–6715.
Ding & Zhou (2020) Ding, X. & Zhou, Z. (2020). Estimation and inference for precision matrices of nonstationary time series. Ann. Statist. 48, 2455–2477.
Ding & Zhou (2021) Ding, X. & Zhou, Z. (2021). Simultaneous sieve inference for time-inhomogeneous nonlinear time series regression. arXiv preprint arXiv:2112.08545 .
Ding & Zhou (2023) Ding, X. & Zhou, Z. (2023). Auto-regressive approximations to non-stationary time series, with inference and applications. The Annals of Statistics 51, 1207 – 1231.
Ding & Zhou (2024) Ding, X. & Zhou, Z. (2024). Supplement to ”On the partial autocorrelation function for locally stationary time series: characterization, estimation and inference” .
Dégerine & Lambert-Lacroix (2003) Dégerine, S. & Lambert-Lacroix, S. (2003). Characterization of the partial autocorrelation function of nonstationary time series. Journal of Multivariate Analysis 87, 46–59.
Goerg (2012) Goerg, G. M. (2012). Testing for white noise against locally stationary alternatives. Statistical Analysis and Data Mining: The ASA Data Science Journal 5, 478–492.
Karmakar et al. (2022) Karmakar, S., Richter, S. & Wu, W. B. (2022). Simultaneous inference for time-varying models. Journal of Econometrics 227, 408–428.
Killick et al. (2020) Killick, R., Knight, M. I., Nason, G. P. & Eckley, I. A. (2020). The local partial autocorrelation function and some applications. Electronic Journal of Statistics 14, 3268 – 3314.
Kley et al. (2019) Kley, T., Preuß, P. & Fryzlewicz, P. (2019). Predictive, finite-sample model choice for time series under stationarity and non-stationarity. Electronic Journal of Statistics 13, 3710–3774.
Ljung & Box (1978) Ljung, G. M. & Box, G. E. (1978). On a measure of lack of fit in time series models. Biometrika 65, 297–303.
Mayer et al. (2020) Mayer, U., Zähle, H. & Zhou, Z. (2020). Functional weak limit theorem for a local empirical process of non-stationary time series and its application. Bernoulli 26, 1891–1911.
Politis et al. (1999) Politis, D., Wolf, D., Romano, J., Wolf, M., Bickel, P., Diggle, P. & Fienberg, S. (1999). Subsampling. Springer Series in Statistics. Springer New York.
Ramsey (1974) Ramsey, F. L. (1974). Characterization of the partial autocorrelation function. The Annals of Statistics , 1296–1301.
Roueff & Sanchez-Perez (2018) Roueff, F. & Sanchez-Perez, A. (2018). Prediction of weakly locally stationary processes by auto-regression. Lat. Am. J. Probab. Math. Stat. 15, 1215–1239.
Shumway & Stoffer (2017) Shumway, R. & Stoffer, D. (2017). Time Series Analysis and Its Applications: With R Examples. Springer Texts in Statistics. Springer International Publishing, 4th ed.
Vogt (2012) Vogt, M. (2012). Nonparametric regression for locally stationary time series. Ann. Statist. 40, 2601–2633.
Wu (2005) Wu, W. (2005). Nonlinear system theory: Another look at dependence. Proc Natl Acad Sci U S A. 40, 14150–14151.
Yuan (2010) Yuan, M. (2010). High dimensional inverse covariance matrix estimation via linear programming. J. Mach. Learn. Res. 11, 2261–2286.
Zhou (2013a) Zhou, Z. (2013a). Heteroscedasticity and autocorrelation robust structural change detection. J. Am. Stat. Assoc. 108, 726–740.
Zhou (2013b) Zhou, Z. (2013b). Inference for non-stationary time-series autoregression. Journal of Time Series Analysis 34, 508–516.
Zhou & Wu (2009) Zhou, Z. & Wu, W. (2009). Local linear quantile estimation for non-stationary time series. Ann. Stat. 37, 2696–2729.
Zhou & Wu (2010) Zhou, Z. & Wu, W. (2010). Simultaneous inference of linear models with time varying coefficents. J.R. Statist. Soc. B 72, 513–531.

\printhistory