Optimal break tests for large linear time series models

Optimal Break Tests for Large Linear Time Series Models

Introduction

This paper uses a growing-dimension asymptotics approach to develop a class of asymptotically valid tests for the stability of infinite or growing dimensional linear regression models, which we will call `large linear models', when the change point is unknown. These tests are motivated by optimality ideas explored in the parametric literature and are indeed shown to possess analogous optimal properties under certain conditions. Our terminology of infinite-dimensional models includes both autoregressive models of infinite-order and nonparametric regression. This paper is thus a contribution to the recently growing interest in making inference on a parameter vector whose dimension is not finite.

Growing-dimension asymptotics have been studied in the statistical literature since at least the important work of huber1973robust, but the sequences of experiments in LeCam1960 constitute an even earlier reference. Our growing-dimensional approach is general enough to include many important models as special cases. For example, we cover the autoregressive model of infinite-order, AR($\infty$), one of the most important models in time series analysis. Since a nonparametric regression function can be represented as an infinite-order linear regression in a linear sieve space, our methods cover this case as well. We allow for a triangular array structure and can thus also accommodate linear regression with growing dimension, but being explicit about the approximation error allows the coverage of true infinite dimensional models such as AR($\infty$) and nonparametric regression. Our focus on linear regression facilitates clarity of exposition, although extensions are feasible to more general settings like nonlinear regression, maximum likelihood, or generalized method of moments (GMM).

While we focus on testing for the presence of a structural break in the AR($\infty$) model, nonparametric regression and growing dimension linear regression, we note that there is a thriving literature on the estimation of high-dimensional change point models, often with LASSO-type methods, e.g. Chan2014; Li2016; lee2016lasso; wang2018high; lee2018oracle; Safikhani2022; Wang2024. This topic has also attracted much attention also in the fixed-dimensional parametric setting, e.g. bai1998estimating, among many others. Our framework is also different from another closely related literature, which explores testing of stability of means or variance-covariances of multiple time series, e.g. aue2009break, cho2015multiple, wang2022inference, li2024detection. In a parametric setting, Andrews1993 considers testing under a GMM framework, which has been extended to the case of a break occurring anywhere in the sample by hidalgo2013testing, and to the systems of equations setting by qu2007estimating. Andrews1994a and elliott2006efficient explore certain optimality properties which form the basis for the tests considered in this paper.

More generally, testing for structural stability is a heavily studied subject and a huge literature has evolved, as reviewed by perron2006dealing,aue2013structural. In the infinite and growing dimensional-setting, however, the literature is not as rich. Wu1993 consider testing for jumps in a regression function using kernels and under a fixed design, vogt2015testing studied break testing in nonparametric kernel regression, and Mohr2019 consider testing in a time series regression model using marked empirical processes of residuals, while Hidalgo2016 use a non-smoothing approach to test the stability of a nonparametric regression function.

To explain our approach further, we view our methodology as approximating an infinite-dimensional testing problem with a sequence of finite-dimensional problems. In case of the nonparametric regression or AR($\infty$), this is based on sieve or series estimation which entails approximating the infinite dimensional primitive with a linear combination of increasingly many basis functions. On the other hand, if the model of interest is an increasing dimensional one, the approach is interpreted as inference on a sequence of increasingly complex models. Because the sequence of basis functions (nonparametric regression), lags (AR($\infty$)) or covariates (increasing dimension regression) that underlies our approach is not of fixed dimension, one cannot simple apply the theory of finite-dimensional parametric problems. Instead, the `usual' parametric test statistics must be adjusted to take into account the increasing dimensional nature of the problem, as shown by Portnoy1988. Murphy1993 used a similar approach in a proportional hazards model while Hong1995 propose an appropriately standardized test statistic that provides a consistent specification test based on series estimation. Fan2001 provides a Wilks Phenonemon type result based on a standardized statistic.

Our testing problem also features a loss of identification under the null of no structural break. Thus, we consider a stochastic process of recentered and renormalized quadratic forms indexed by the unidentified parameter that captures the change point as a fraction of the sample size. We establish the weak convergence of this stochastic process as both the sample size $T$ and the number $p $ of restrictions grow to infinity simultaneously and construct test statistics based on a weighted exponential transformation motivated by the optimal tests in Andrews1994a. Under an i.i.d. Gaussian design, we establish a weighted power optimality for our class of tests in the same spirit.

Optimality typically requires a finite number of directions in which the null may be violated, see e.g. Section 14.6 of Lehmann2006. To account for the infinite-dimensional aspect of the problem, we adopt the spirit of sieve regression and assume that we can direct our test to have reasonable power against an increasing number of approximating directions, while the remaining `tail' directions can be neglected. Among the approximating directions, we assign weights inspired by Andrews1994a. In parametric structural break testing, elliott2006efficient found that most optimal tests are similar as long as they are optimal against alternatives of similar average deviation from the null. When $p$ grows to infinity, small changes in each coefficient can aggregate, so that individually negligible changes can add up to a non-negligible change as a whole. Indeed, we show that the determinant of our local power is the average of all the changes.

In the general model without the i.i.d. Gaussianity assumption, we observe that an important feature of the limiting process is the presence of a nonlinear high order long run variance (HLV) involving the errors and regressors in a nonlinear fashion, Gupta2023. We also show that the sequential limit fails to uncover this nonlinear serial correlation, where the limits are taken sequentially first as $T\rightarrow \infty $ and then as $p\rightarrow \infty .$ We robustify our test against this HLV factor using a random scaling method as opposed to resorting to consistent estimation. This yields a pivotal test and is in the spirit of the self-normalization approach, see e.g. shao2010testing, Zhang2018. In addition, we incorporate a bootstrap bias correction that controls size effectively in finite samples while preserving power.

We also examine the finite-sample performance of our test in a simulation study. We find that our test exhibits very good size control and high power against a variety of alternatives, both in nonparametric regression and long AR estimation. Importantly, we find that failure to correct for the nonlinear serial correlation that we have discovered can lead to seriously distorted sizes. We show that this distortion can be particularly severe as $p,T$ increase to infinity, and the correction we propose dramatically improves matters in many cases. In an empirical illustration with real data we examine the relationship between oil prices and economic output and find that our corrections can yield different conclusions to existing methods.

The next section introduces the class of optimal tests that motivate our approach, sets out the mise en sc\`ene of our structural break model and establishes basic asymptotic results and the optimality property of our tests. In Section (ref) we show how to construct HLV-robust and bias corrected versions of our tests and establish their asymptotic properties, while Section (ref) shows that taking a sequential limit approach would hide the HLV factor that our simultaneous limit theory reveals. In Section (ref) we demonstrate the desirable finite sample properties of our tests in a set of Monte Carlo simulations. Finally, we illustrate our tests on real data in Section (ref).

Optimal Testing

We develop an optimal test for the presence of a structural break in infinite-order regression such as the conditional mean of nonparametric regression models with iid Gaussian errors in terms of a weighted average power. The growing dimensional regression model is a special case of our model. It generalizes conventional optimal structural break testing, which Andrews1994a proposed for finite-dimensional regression. Specifically, we first characterize the optimal test under:

which defines a sequence of models where $\psi _{Tt}=\psi -1\left\{ t/T\leq \psi \right\} $ for $\psi \in \Psi=\left[c_1,c_2\right] \subset (0,1)$, $\{\nu _{t}\}$ is an i.i.d. sequence of standard normals and is independent of the time series of infinite-dimensional vectors $\left\{ \boldsymbol{x}_{t}\right\} $, while $p$ grows with the sample size $T$. We drop the subscript $T$ to ease notation hereafter unless it is necessary to avoid confusion. The null hypothesis of interest is given by $$H_0 : (b_{1},b_{2},...)'=0, \quad \text{for any } \psi \in \Psi. $$

For given weight functions for $\psi$ and $\boldsymbol{b}= (b_{1},b_{2},...)'$, which we denote by $J\left( \psi \right) $ and $Q\left( \boldsymbol{b}\right) $ respectively, the weighted average power of a test $\mathscr{T}_{T}$ over this class of local alternative models is written as

where $Y$, $ \boldsymbol{X}$, and $p_{\boldsymbol{b},\psi }\left(Y|\boldsymbol{X}\right)$ denote the collections of $\left\{ y_{t}\right\} $, $\left\{ \boldsymbol{x}_{t}^{\prime }\right\} ,$ and the conditional density, respectively, and the equality holds by Fubini's theorem. The final expression reveals that an optimal test can be constructed so as to maximize the weighted average likelihood $\int \int p_{\boldsymbol{b},\psi }\left( Y|\boldsymbol{X}\right) dQ\left( \boldsymbol{b}\right) dJ\left( \psi \right) $ for any given probability measures $Q$ and $J$. In this alternative representation, the alternative hypothesis is a simple hypothesis. That is, the alternative hypotheses can be written as \[H_{A}:\int \int \phi \left( \boldsymbol{x}_{t}^{\prime }\gamma _{0}+\left( T\sqrt{p}\right) ^{-1/2}\psi _{t}\boldsymbol{x}_{t}^{\prime }\boldsymbol{b},\sigma_{0}\right) dQ\left( \boldsymbol{b}\right) dJ\left( \psi \right) , \] where $\phi \left( \mu ,\sigma \right) $ denotes the normal density function with mean $\mu $ and variance $\sigma ^{2}$. Assuming $\gamma _{0}=0$ (without loss of generality) yields $y_{t}=\sigma _{0}\nu _{t}$ under the null. Then, the likelihood ratio (LR) statistic is given by

As noted by Andrews1994a, this is a useful feature since optimality can be established for the likelihood ratio test by virtue of the Neyman-Pearson lemma.

Thus, it is important to specify the weights on $\boldsymbol{b}$ and $ \psi $ in an appropriate fashion. A particular choice of the weighting function was proposed by Wald1943 for a fixed $\psi $ and extended by Andrews1994a and Song2009 for the case where $\psi $ is not identified. The latter papers consider a finite number of restrictions on functions of $\psi $ and the weight function assigns constant weight on the ellipses in the parameter space of $ \boldsymbol{b} $ that are equally difficult to detect for each $\psi$. As noted by Andrews1994a, every locally most powerful invariant test has to exhibit a constant power over these ellipses under the parameterization in (ref). Our testing problem concerns an infinite-dimensional restriction that is also a function of $\psi$.

\sloppy Building on the observation by Lehmann2006 (Section 14.6) that one cannot obtain reasonable local power across all families of distributions when the restriction concerns an infinite-dimensional parameter, we consider weighting functions that coincide with these ellipses only up to a truncation parameter, which is set as $p$ in (ref). Specifically, let $x_{t}$ stand for the first $p$-elements of $\boldsymbol{x}_{t}$ and let $w_{t}$ denote the remaining elements so that $ \boldsymbol{x}_{t}=(x_{t}^{\prime },w_{t}^{\prime })^{\prime }$ and partition $\boldsymbol{b}=\left( \boldsymbol{b}_{1}^{\prime },\boldsymbol{b} _{2}^{\prime }\right) ^{\prime }$ accordingly. Recall that we drop $T$ in the notation. Also, write $ Q\left( \boldsymbol{b}\right) =Q_{1}\left( \boldsymbol{b}_{1}\right) Q_{2}\left( \boldsymbol{b}_{2}|\boldsymbol{b}_{1}\right) , $ where $Q_{1}\left( \boldsymbol{b}_{1}\right) $ is the $p$-dimensional centered multivariate normal distribution function with covariance matrix $ \mathfrak{c}\sigma_0^{2}\left( \psi \left( 1-\psi \right) E{x}_{t}{x}_{t}^{\prime }\right) ^{-1}$ for some finite $\mathfrak{c}$ and $Q_{2}\left( \boldsymbol{b}_{2}| \boldsymbol{b}_{1}\right) $ is a conditional distribution over the support $ \left\{ \boldsymbol{b}_{2}:\left\Vert \boldsymbol{b}_{2}\right\Vert _{2}\leq K\ \text{and }\sqrt{p}\left\Vert \boldsymbol{b}_{2}\right\Vert _{2}=o\left( 1\right) \right\}, $ for $\boldsymbol{b}_{2}$ given $\boldsymbol{b}_{1}$ and some finite $K$. The growing dimensional model sets $\boldsymbol{b}_2 = 0$. The weight $J$ for $\psi $ may be set as uniform on the interval $\Psi $ unless there is a specific set of break points of interest.

The restriction on the decay rate of the infinite-dimensional parameter $\boldsymbol{b}_{2}$ arises naturally from the features of the particular model under consideration. For instance, berk1974consistent and Kuersteiner2005 constrain the autoregressive coefficients to meet the decay rate of $\sum_{j=1}^{\infty }j\left\vert \gamma _{j}\right\vert <\infty $ for the stationarity and moment conditions of the process and consistent estimation of the spectral density at the origin. On the other hand, in the nonparametric regression context, the smoothness condition on the regression functions can be imposed by constraining the decay rate of the sieve coefficients.

To illustrate, consider the Fourier representation $f\left( x\right) =\sum_{j=1}^{\infty }\gamma _{j}\exp \left( ijx\right) $, for which the $m$-th derivative is given by $f^{\left( m\right) }\left( x\right) =\sum_{j=1}^{\infty }\gamma _{j}\left( ij\right) ^{m}\exp \left( ijx\right). $ Thus, it is natural to assume that $ \sum_{j=1}^{\infty }\left\vert \gamma _{j}\right\vert j^{m}<\infty $ for $ f^{\left( m\right) }\left( X\right) $ to be well-defined a.s. This also implies that $\gamma _{j}=o\left( j^{-m-1}\right) $. These observations are reflected in our construction of the weight function. The preceding specification of local models with truncation $p$ enables us to construct a class of computationally tractable tests that are asymptotically equivalent to the LR test. We develop such an asymptotically optimal test in the following subsection, extending the finite-dimensional parametric optimality results of Andrews1994a.

Structural Break in Growing and Infinite-Dimensional Regression

We consider the linear regression model of growing dimension with a potential structural break at $ t=[T \psi] $,

where $\psi\in\Psi$ is unknown, $x_t$ is a $p$-dimensional observable vector, $\kappa_1, \kappa_2$ are $p$-dimensional unknown parameters, and $p\rightarrow \infty $ as $T\rightarrow \infty $. This is the truncated model we employ to construct our test statistic and $\kappa_2$ corresponds to $\boldsymbol{b}_1$ so that we construct a test for the null hypothesis \[ \mathcal{H}_0 : \kappa_2 =0 \quad \text{for any } \psi \in \Psi. \]

We will allow for general increasing dimension and nonparametric time series regressions, $E\left( y_{t}| \mathfrak{G}_{t-1}\right) $, where $\mathfrak{G}_{t-1}$ denotes the filtration up to time $t-1$, thus generalizing the Gaussian model (ref). Let $\nu _{t}=y_{t}-E\left( y_{t}|\mathfrak{G}_{t-1}\right) $, $\gamma_T$ be the best linear predictor of $y_{t}$ given $x_{t}$, and $r_{Tt}=E\left( y_{t}|\mathfrak{G}_{t-1}\right) -x_{Tt}^{\prime }\gamma_T $. Then, $u_{t}=r_{Tt}+\nu _{t}$, i.e. the error term comprises of the regression error in addition to the approximation error. Here we use the subscript $T$ to emphasize the the array structure temporarily.

For example, for growing dimensional regression the model features an increasing number of observable covariates $x_{t}$ as the sample size diverges, a modelling approach that dates back to the pioneering work of huber1973robust. In the case of infinite order autoregression, we include the lagged dependent variables, $\left\{y_{t-j}\right\}_{j\geq 1}$, and our approach indicates a truncation at $p$ lags, $p\rightarrow\infty$, as suggested in the classical contributions of berk1974consistent and Shibata1980, among others. Finally, in the case of nonparametric series regression, the model features a finite number of observable covariates $z_{t}$ and $ x_{Tt}:=x_{T}\left( z_{t}\right) :{\mathbb{R}}^{k}\mapsto {\mathbb{R}}^{p}$ are transformations via known basis functions.

Let $x_{t}\left( \psi \right) :=\left(x_{t}^{\prime },x_{t}^{\prime }1\left\{ t/T>\psi \right\} \right) ^{\prime}$ and denote by $ \hat{\kappa }\left( \psi \right) $ and $ \hat{u}_t\left( \psi \right) $ the OLS estimate and the OLS residuals, respectively, from the linear regression with a fixed $\psi$. Define $ \hat{P}\left( \psi \right) ={T}^{-1}\sum_{t=1}^{T}x_{t}\left( \psi \right) x_{t}\left( \psi \right) ^{\prime }, $ and let $\hat{\Xi}\left( \psi \right) $ denote an estimator of $ E\nu _{t}^{2}x_{t}\left( \psi \right) x_{t}^{\prime }\left( \psi \right) $. One can use $T^{-1}\sum_{t=1}^{T}x_{t}\left( \psi \right) x_{t}\left( \psi \right) ^{\prime }\hat{u}_{t}\left( \psi \right) ^{2}$ (the Eicker-White formula) or even $\hat{\sigma}\left( \psi \right) ^{2} \hat{P}\left( \psi \right) $, where $\hat{\sigma}^{2}\left( \psi \right) =T^{-1}\sum_{t=1}^{T}\hat{u}_{t}\left( \psi \right) ^{2}$, if conditional homoskedasticity obtains. Define the heteroskedasticity robust Wald process as

where $S=\left( 0_{p\times p}:I_{p}\right) $, with $I_p$ the $p\times p$ identity matrix. To account for the growth of $ p $ as a function of $ T $, we center and rescale $W_T(\psi)$ as

The idea is that this standardization yields a proper limit as $p$ and $T$ diverge to infinity simultaneously, motivated by the intuition that a correctly centered and scaled sequence of chi-squared random variables converges to a normal variate in distribution. It is motivated by Portnoy1988, DeJong1994, Hong1995, Fan2001, Gupta2023, for example. They established a CLT for the centered and rescaled statistic $\mathcal{Z}_{T}\left (\psi \right )$ for the case where $\psi$ is fixed at a constant, under various sampling conditions.

We extend the CLT for time series data in Gupta2023 to a functional CLT (FCLT) by establishing the weak convergence of ${\mathcal{Z}}_{T}\left (\psi \right )$ as a process indexed by $\psi\in \Psi$ in order to construct an optimal test for the structural break with an unknown break date. Define the process $$ \mathcal{Z}(\psi )={\psi^{-1}G(\psi )+(1-\psi)^{-1}\bar { G}\left (\psi \right )-G(1)},$$ where $\left (G\left (\psi \right ),\bar {G}\left (\psi \right )\right )^{\prime }$ is a bivariate Gaussian process with covariance kernel

Observe that the marginal distribution of $ \mathcal{Z}(\psi ) $ is standard normal.

\sloppy Under the presence of strong high-order serial dependence, the limit of the finite-dimensional distribution of $ \mathcal{Z}_T (\psi)$ features the HLV factor as shown by Gupta2023. For convenience, we let $\varsigma_{t}=\Xi^{-1/2}x_{t}\nu_{t}$, where $\Xi$ is defined in Assumption (ref), $\delta_t(\psi)=\left({1}\left\{t/n\leq \psi\right\}-\psi\right)/\sqrt{\psi(1-\psi)}$ and define the HLV factor as

where $q_t(\psi)=(Tp)^{-1/2}\delta_t(\psi)\varsigma_t'\sum_{s<t}\delta_s(\psi)\varsigma_s$, and any limit stated as `$T\rightarrow\infty$' is taken as both $ T $ and $ p $ grow to infinity simultaneously. A sequential limit is considered in Section (ref). Since the variable $q_t(\psi)$ is a nonlinear transformation of the underlying primitive random variables $x_t$ and $\nu_t$ and their entire past, its variance depends on $t$ and contains higher-order autocovariances of the primitive variables or $\varsigma_t$. Higher-order dependence of a martingale difference sequence (mds) process is a common feature of many time series as in e.g. the ARCH process. In this sense it differs from the familiar long-run variance encountered in fixed dimensional time series regressions, which is the limiting sum of all the autocovariances.

Throughout let $C$ denote a generic finite constant, independent of $T$.

Introduce the $p\times p$ non-random matrix sequences $P$ and $ \Xi $ that satisfy Assumption (ref) below and let

Define $\left\Vert A\right\Vert =\left\{ \overline{\text{eig} }(A^{\prime }A)\right\} ^{\frac{1}{2}}$ for a generic matrix $A$, where $\underline{ \text{eig} }$ (respectively $\overline{\text{eig} }$) denotes the smallest (largest) eigenvalue of a symmetric nonnegative definite matrix.

Primitive conditions and expressions for $\tau _{p}$ and $\zeta_p$ are given in Propositions B1 and B2 of Gupta2023. Recall that the eigenvalues of the Kronecker product of two symmetric matrices are the products of their eigenvalues, and $ \psi $ is bounded away from zero and one. Thus, $P(\psi)$ and $\Xi(\psi)$ inherit the eigenvalue restrictions on $P$ and $\Xi$ in Assumption (ref) $(ii)$ and $(iii)$, up to positive constants.

Let $ \mathfrak{H}_t $ denote a filtration for $ \varsigma_{t} $, $\Phi_t=E\left(\varsigma_{t} \varsigma_{t}'|\mathfrak{H}_{t-1}\right)$, and $\Theta_s=\sum_{t_1=1}^{s-1}\sum_{t_2=1}^{s-1} \varsigma_{t_1} \varsigma_{t_2}'$. The filtration $ \mathfrak{H}_t $ need not be $ \mathfrak{G}_t $ but a simpler one as long as it makes $ \varsigma_t $ an mds. Indeed, some conditions may be easier to verify under simpler filtrations. The next assumption introduces the HLV factor $\upomega$ formally.

As we discussed subsequently to (ref), this quantity in the asymptotic distribution differs from the conventional long-run variance and consists of fourth-order cross-moments of a time series. Then, a resampling algorithm matching only up to second order moments such as the wild bootstrap would not be valid. Also note that if the $v_t$ are i.i.d. Gaussian then $\upomega=1$.

For mean zero random variables $a_{1i},a_{2j},a_{3k},a_{4l}$, let $\mathrm{cum}_{ijkl}\left(a_{1i},a_{2j},a_{3k},a_{4l}\right)$ denote the fourth cumulant.

This assumption controls the temporal dependence in $\left\{x_t\nu_t\right\}$ and is discussed in Andrews1991a, for example, wherein sufficient conditions for it to hold are also provided. The next two theorems establish the stochastic equicontinuity of the process $\mathcal{Z}_T(\psi)$ and its weak convergence. Let `$\Rightarrow$' denote weak convergence in $\ell^{\infty}(\Psi)$.

For local power properties, let the sequence of local alternatives be

where $\varrho $ is a unit length $p\times 1$ vector. Due to this unit length condition, the power of $p$ becomes $1/4$ from $-1/4$ in (ref). This sequence converges to the null slower than the usual $1/\sqrt{T}$ parametric rate due to the inflation caused by the $p^{1/4}$ factor. This slower rate is the price to pay for a nonparametric approach.

Asymptotic Optimality

To construct our test statistics, we exploit an approach to testing that was initiated by Andrews1994a for parametric models and has been widely applied e.g. Song2009. This involves constructing test statistics by weighted exponential transformations of an underlying process. These are motivated by the fact that, under certain simplifying assumptions, tests constructed in this fashion have an optimality property, as we explore in Section (ref). In the parametric setting, elliott2006efficient showed that many of the seemingly different optimal tests are in fact asymptotically equivalent within a large class of local break models and what matters is the weighted average size of the breaks, regardless of the number of breaks. These references do not treat the stability of nonparametric or growing dimensional models.

Accordingly, define the following class of test statistics:

for a positive $\mathfrak{c}$ and a bounded weight function $J(\cdot)$ such that $\int_{\Psi }dJ\left( \psi \right) =1$. These are Andrews1994a-style test statistics but correctly constructed to account for $p\rightarrow\infty$ as well as the centered and rescaled Wald process $\mathcal{Z}_{T} \left( \psi \right)$ that we introduced above.

The optimality of the $Exp{\mathcal{Z}}_{T}$ test under the Gaussian design is formally given by the following equivalence theorem. We will elaborate on the robustification against the HLV factor in the subsequent section.

To shed light on our optimality, we elaborate the relation of our $Exp{\mathcal{Z}}_{T}$ test statistic to the optimal test statistic Exp-W introduced by Andrews1994a for the parametric model, viz. $\left( 1+ \mathfrak{c}\right) ^{-\frac{p}{2}}\int_{\Psi }\func{exp}\left( \frac{1}{2} {\mathfrak{c}}/({1+\mathfrak{c})}W_{T}\left( \psi \right) \right) dJ\left( \psi \right) $, where $W_{T}\left( \psi \right) $ is the Wald statistic defined in ((ref)). Specifically, (ref) in the proof of Theorem (ref) shows that

under both the null and local alternatives. For a moderate value of $ \mathfrak{c}$, our $Exp{\mathcal{Z}}_{T}$ statistic can be approximated by an exponential statistic Exp-W with some $\mathfrak{c}^{\prime }$ by a change-of-variables. On the other hand, the conventional supremum statistic for the structural break testing can also be represented by a limit case of our $Exp{\mathcal{Z}}_{T}$ since $c^{-1}\log \sum_{i=1}^n exp(cx_i) \to \max_i x_i $ as $c\to \infty$. This is a different feature from Andrews1994a.

A Class of Optimality Motivated HLV-Robust and Pivotal Tests

HLV Robustification of $ \mathcal{Z}_T(\psi) $

As observed above, the HLV term $\upomega$ is a type of long run variance and hence can handled with techniques for this kind of problem. Heteroskedasticity and autocorrelation consistent/robust (HAC/HAR) inference or the estimation of the spectral density at frequency zero has a long history in time series literature, see e.g. priestley1988spectral,newey1987simple,Andrews1991a,kiefer2002heteroskedasticity,shao2015self. Noting that the asymptotic covariance kernel of the process $ \mathcal{Z}_T(\psi) $ depends on HLV, we utilize the fixed-bandwidth kernel approach to obtain an asymptotically pivotal test. For this purpose, $ k(\cdot) $ denotes a kernel function that satisfies Assumption (ref).

Since $\nu_{t}$ and $\Xi$ are not directly observable in practice, we replace them with the least squares-based estimates as in Section (ref) and introduce $\ell _t = \left( np\right)^{-1/2} x_{t}'\hat{\Xi}^{-1} \hat{\nu}_{t} \text{\ensuremath{\sum_{s=1}^{t-1}x_{s}\hat{\nu}_{s}}} $ and its demeaned version, $\bar{\ell}_{t}= \ell_t - T^{-1}\sum_{t=2}^{T} \ell_t$. A feasible kernel weighted estimate of $\upomega$ is then given by

Gupta2023 observed that a fixed bandwidth approach, as in Sun2014, works well when dealing with corrections for long-run variances of nonlinear transformations of primitive variables, which is what $\upomega$ captures. Our estimator is accordingly based on the weighting function $ K_{h}\left(r,s\right) = k\left(h\left(r-s\right)\right)$, where $ h=1/b $ and $k\left(u\right)=\left(1-\left|u\right|\right)^{h}1\left\{ \left|u\right|<1\right\} $, using the Bartlett kernel with $h=1$. This is what is termed by Sun2014 as the sharp kernel estimator. In the theorems presented below, we characterize the joint weak limit of $ \hat{\upomega} $ and $ \mathcal{Z}_T (\psi) $, but also incorporate a bootstrap bias correction.

Boostrap Bias Correction and Pivotalization of $ \mathcal{Z}_T(\psi) $

Observe that the degrees of freedom $p$ provide a first-order asymptotically correct centering for $ W_T (\psi) $ but Gupta2023 show that this might induce bias in finite samples and propose a bootstrap bias correction. This is done via the null-imposed wild bootstrap by generating

where $ \upsilon_t $ is an iid sequence of centered and normalized variables, e.g. Rademacher variables, to compute $ \mathcal{Z}_T^{\star} (\psi) $. The imposition of the null helps enhance the power of the test. Repeat $ B $ times to obtain $ \bar{\mathcal{Z}}_T^{\star} (\psi) = B^{-1} \sum_j^{B} \mathcal{Z}_T^{\star, j} (\psi) $, the bootstrap estimate of the bias. It is worth commenting that this bootstrap may not properly approximate the distribution of the test statistic even asymptotically, except the first moment that we utilize in our correction. This is due to the dependence of the asymptotic variance of the test statistic on the higher order moments of the sampling distribution.

Then, define the bias corrected and HLV-robustified version of $\mathcal{Z}_T (\psi)$:

Now, let the superscript $ \star $ indicate the bootstrap analogue and $\mathscr{K}_{h}\left(r,s\right) =K_{h}\left(r,s\right)-\int_{0}^{1}K_{h}\left(\varrho,s\right)d\varrho-\int_{0}^{1}K_{h}\left(r,\varrho\right)d\varrho +\int_{0}^{1}\int_{0}^{1}K_{h}\left(\varrho_{1},\varrho_{2}\right)d\varrho_{1}d\varrho_{2} $. Following a standard type of assumption on the kernel function we establishe pivotality of $\mathscr{H}_T^b (\psi)$.

We introduce a class of HLV robust and pivotalized weighted exponential statistics based on our approach that encompasses nonparametric and growing dimensional models. Define

again for a positive $\mathfrak{c}$ and a bounded weight function $J(\cdot)$ such that $\int_{\Psi }dJ\left( \psi \right) =1$, and for which we set $ \hat{\kappa}_1 = \hat{\kappa}_1 (\hat{\psi}) $ and $ \hat{\nu}_t =\hat{\nu}_t (\hat{\psi})$ in (ref), where $ \hat{\psi} = \arg \min \hat{\sigma}(\psi) $ and $ \hat{\sigma}(\psi) $ is the sum of squared residuals as defined in Section (ref). Since

we may extend the definition of $ Exp{\mathscr{H}}_T (\mathfrak{c} )$ for $\mathfrak{c} \in [0,+\infty]$.

The asymptotic distributions under both hypotheses follow from the continuous mapping theorem given the weak convergences of the stochastic process $\mathcal{Z}_T(\psi)$ on $ \Psi$ that are established in Theorems (ref) and (ref). The main result is stated in the following theorem.

The noncentrality term is also positive for any $ \psi \in \Psi $ to make the test nontrivial.

Sequential Limit

This section elaborates a sequential limit where first $ T \to \infty $ and then $ p \to \infty $, and the relation between the fixed $ p $ critical values and our critical values for the two extreme limiting cases of the weighted exponential tests. We show that the HLV factor $\upomega$ is wiped out in the sequential limit. First recall that, for each $ p $, the limit has been obtained by Andrews1993. Indeed, fixing $p$, Andrews1993 has showed that the $W_T(\psi)={\mathcal{L}_{T}}(\psi )+o_p(1)$ uniformly in $\psi\in\Psi$, and the weak limit of the process

is

where $B_{p}$ stands for the $p$-dimensional standard Brownian motion and thus $\mathcal{W}_{p}$ is the standardized tied-down Bessel process of degree $p$. For each $\psi $, $\mathcal{W}_{p}\left( \psi \right) $ is distributed as a Chi-square with $p$ degrees of freedom, thus with mean $ p$ and variance $2p$. This is also a pivotal process. The following proposition derives the second limit in the sequential limit, i.e. as $ p \to \infty $, showing that it hides the HLV factor $\upomega$.

A consequence of Proposition (ref) is that Andrews' (1993) critical values, say $ {\alpha }$ are valid after the transformation

due to Lemma 21.2 in VanderVaart1998, assuming that the limit distribution function of our test statistic is continuous at the $\alpha $ -level critical value. Let $c_{\alpha }^{\ast }$ and $c_{\alpha }^{p}$ be the solutions of $ \Pr \left\{ \sup_{\psi }\mathcal{Z}\left( \psi \right) >c\right\} =\alpha \ \ \text{and } \ \Pr \left\{ \sup_{\psi }\mathcal{W}_{p}\left( \psi \right) >c\right\} =\alpha , $ respectively. Then, as

as $p\rightarrow \infty ,$ by Proposition (ref), we conclude from Lemma 21.2 in VanderVaart1998 that

The same argument holds true for the average test. For general weighted exponential tests, the transformation for the critical values is more involved than in (ref) as given by the relationship (ref) below with some approximation error.

Monte Carlo

In this section we examine the size and power performance of $Exp{\mathscr{H}}_{T}(\mathfrak{c})$ for the presence of a structural break at an unknown date, focusing on a special case of the exp test, viz. $Exp{\mathscr{H}}_{T} := \lim_{\mathfrak{c}\rightarrow \infty}Exp{\mathscr{H}}_T(\mathfrak{c}) $, see ((ref)). Recall that it is equivalent to $\sup_{\psi\in\Psi}\mathscr{H}_{T}^{b}\left(\psi\right)$. We compare it with $Exp{\mathcal{Z}}_{T} := \lim_{\mathfrak{c}\rightarrow \infty}Exp{\mathcal{Z}}_T(\mathfrak{c}) $ and $supW_T :=\sup_{\psi\in\Psi}W_{T}\left(\psi\right)$, where the critical values for the last test come from Andrews1993. The experiments for size cover $T=250,500,$ while for power we present results for $T=400$. Finally, we set $\Psi=[0.15,0.85]$.

We consider two examples: multiple regression and infinite order autoregression or the ARMA model, for which we generate the error from a bounded ARCH process

where $ f(x) = x^2 1\{|x|\leq b\} +b^2 1\{|x|> b\} $, $\bar\eta_{t}=\left(\eta_{t}-E\eta_{t}\right)/\sqrt{var\left(\eta_{t}\right)}$, and $\left\{ \eta_{t}\right\} $ is an iid sequence from the normal mixture distributions of type 1 and 2 in Marron1992. We fix $b=2.5$. Their mixture type 1 is the standard normal. For a standard normal vector $\left(Z_{1},...,Z_{k}\right)$ and multinomial vector $\left(d_{1},...,d_{k}\right)$ with probability $\left(1/5,1/5,3/5\right)$, the mixture type 2 skewed unimodal variate is $\eta_{t}=Z_{1}d_{1}+\left(2Z_{2}/3+1/2\right)d_{2}+d_{3}\left(5Z_{3}/9+13/12\right)$. We also experimented with their mixture type 3 but do not report this here because the results are similar.

More specifically, for the multiple regression, the regressors $ x_t $ consist of independent AR(1) processes with coefficient $ \alpha_x $ and the ARCH innovation as in (ref) and their lags of order up to 3. That is, we consider the distributed lag (DL) model with growing numbers of variables. The first five elements of the coefficients are set as $d_0\left(5^{-1/2},...,5^{-1/2}\right)p^{1/4}T^{-1/2}$ and the others as zeros. When there is a break, all the values become zero after the break so that the value $ d_0 $ controls the magnitude of the change. We vary $\alpha \in \{0.3,0.55\}$, and $p \in \{5,9,13\},$ to examine the effect of the dimension on our tests. For the infinite order AR regression, we generate the sample from the MA(1) model $ y_{t}=\nu_{t}+\theta1\left\{ t\leq\mu\right\} \nu_{t-1} $, $\mu=T\psi$, and estimate the AR($ p $) model with $p=9$ for $ T = 250 $ and $ p = 18 $ for $ T = 500 $. Note that the coefficient $\theta$ represents the size of the jump as well as the MA coefficient.

The vertical partitions in Table (ref) represent different values of $p$, increasing from left to right. The traditional sup test $supW_T$ over-rejects under both $T=250$ and $500$, with performance getting poorer as $p$ increases. The correction for $p$ immediately improves matters, with $Exp{\mathcal{Z}}_T $ still over-rejecting but to a lesser degree. Our recommended statistic $Exp{\mathscr{H}}_T$ performs best overall, and achieves quite acceptable size performance when the sample size is bigger with $T=500$ even when $p$ is moderate. Much the same lesson is learnt from Table (ref), that is, failure to correct for the HLV factor $\upomega$ or the finite sample bias leads to size distorted tests. We note that the $supW_T$ and $Exp{\mathcal{Z}}_T$ tests are severely undersized while our recommended statistic $Exp{\mathscr{H}}_T$ maintains reasonable rejection rates that are much closer to the nominal $5$ % level. Here, we only report the case of $T=500$ to save space since the rejection rates for $supW_T$ and $Exp{\mathcal{Z}}_T$ are more than 0.5 when the $T=250$.

Finally, each panel of Table (ref) displays power for both regressions, showing a vertical partition that corresponds to increasing break magnitude from left to right. Specifically, the partitions correspond to $d_0=1,5,10$ for the multiple regression and $\theta = 0.4,0.6,0.8$ for ARMA. These parameters control the size of the change of the coefficients after the break. Within each partition the parameters vary lexicographically as $(\eta, \alpha)$. We observe that the power of $Exp{\mathscr{H}}_T$ and the other statistics grow as the magnitude of the structural break increases. Across different tests, there is no dominant test in terms of power. It may reflect large size distortions of the $supW_T$ and $Exp{\mathcal{Z}}_T$ tests, but our recommended statistic $Exp{\mathscr{H}}_T$ does not unduly sacrifice power.

Empirical Application

We apply our test to examine the stability of the dynamics of economic output variables and the oil price–output nexus, originally examined in Hamilton2003, using an autoregressive distributed lag model (ADL) of order $(p,p)$. We employ quarterly U.S. macroeconomic data, where the growth rate of chain-weighted real GDP serves as the measure of real activity. Oil price dynamics are captured through the nominal crude oil producer price index, not seasonally adjusted.

Following Hamilton's formulation, we consider three oil price measures: (i) the quarterly percentage change in oil prices, denoted $o_t$; (ii) a one-sided transformation capturing only positive changes, $o_t^+ = o_t \cdot \mathbf{1}\{o_t > 0\}$; and (iii) the net oil price increase, $o_t^T$, which equals the excess of log oil prices in period $t$ over the peak value observed in the prior 12 months, or zero otherwise. The sample spans from 1949:Q1 to 2019:Q4 and is obtained from the FRED database of the St. Louis Federal Reserve.

As a sensitivity measure, we extend the analysis to monthly data using the industrial production (IP) index as an alternative measure of economic activity. Given the monthly frequency, we employ the ADL$(12,12)$ and ADL$(18,18)$ models, which correspond to one-year and 18-month lag structures, respectively. The shift to monthly data substantially increases the number of restrictions involved in testing for structural change. The number of restrictions vary from 13 in AR(12) and 25 in ADL(12,12), to 19 in AR(18) and 37 in ADL(18,18), while the number is at most 13 in case of the quarterly series. This growth in dimensionality has important implications for inference. Notably, critical values for the standard $supW_T$ test are not tabulated for degrees of freedom exceeding 20 in Andrews1993 or Hansen1997a, further motivating the use of our alternative procedure.

We first test for structural instability in the univariate series using AR$(p)$ models. Subsequently, we estimate ADL$(p,p)$ regressions of GDP growth or industrial production index on oil prices using the three oil price specifications, and examine potential structural breaks in these relationships. In line with Hamilton2003, we do not fix a specific break date but adopt $Exp{\mathscr{H}}_T $ statistics evaluated over a trimmed interval $\Psi = [0.15, 0.85]$ as in Section (ref). And we compare it with $supW_T$, where the latter's critical values are computed using the R function provided in Hansen1997a and they are not accounted for the issues arising from large $p$.

Table (ref) reports the p-values of both tests for various specifications. The results reveal divergent inference: while the standard $supW_T$ test supports the existence of structural breaks in both AR and ADL regressions, the robust $Exp{\mathscr{H}}_T$ test yields mixed evidence, suggesting that conventional methods may overstate the presence of instability in this context. The results for the full sample mirror the patterns observed in the Monte Carlo simulation. However, the discrepancies between the robust and conventional test statistics are magnified in the settings of larger $p$s, reinforcing concerns that standard procedures may yield misleading inferences in high-dimensional applications.

We also conduct subsample analyses to assess the temporal robustness of our findings to account for potential structural shifts in macroeconomic dynamics, such as changes in monetary policy regimes or energy market structures. Specifically, we consider two subsamples: Sub 1 starts after the 1980 oil shock while SS2 ends before the 2007 oil shock. As for the lag order $ p $, we try both $ p=4 $ and $ p=6 $ for GDP dynamice and $p=12$ and $p=18$ for IP dynamics as before. For each subsample, we re-estimate the AR and ADL models and apply both the standard and robust structural break tests.

The results are mostly similar to those observed in the full sample. The standard Wald test supports the presence of instability in the oil-output relations most of the time, while our test tells mixed story depending on the specifications. Also, the instabilities in the dynamics of GDP growth rate and industrial production have become substantially weakened. Contrary to the expectation that structural instability might be more pronounced in the earlier period, the p-values across subsamples do not exhibit a systematic pattern. In several cases, the robust test detects weaker evidence of breaks in SS1 relative to SS2, suggesting that the magnitude and nature of instability may vary by specification and measure of oil prices. These findings underscore the importance of flexible testing procedures that do not presuppose the direction or timing of instability. They also caution against broad characterizations of macroeconomic stability across periods without formal statistical validation.

Conclusion

We have developed a class of optimal tests for structural breaks occurring at unknown dates in infinite and growing-dimensional linear time series models. Extending the exponential-average framework of Andrews1994 to the large-dimensional setting, we established average power optimality under Gaussian designs and derived a functional central limit theorem that remains valid when both the sample size and the number of restrictions diverge. The limiting process features a nonlinear high-order long-run variance term that captures interactions between serial dependence and dimensional growth, a phenomenon absent in fixed-dimensional models.

We proposed a random-scaling correction and a bootstrap bias adjustment that jointly yield a pivotal, HLV-robust test. The resulting statistic maintains correct asymptotic size and performs well in finite samples. Monte Carlo experiments demonstrate substantial gains in size control relative to conventional supremum-type and exponential tests, particularly when the number of regressors is large. The framework integrates optimal testing principles with growing-dimension asymptotics and provides a basis for more robust and interpretable inference than standard tests, highlighting its potential value in high-dimensional data analysis where conventional asymptotics are unreliable.