Detecting multiple change points in linear models with heteroscedastic errors

Detecting multiple change points in linear models with heteroscedasticity

Introduction

Linear models are widely used for causal inference and out-of-sample prediction problems with time series data, including macroeconomic forecasting, asset pricing, and portfolio optimization. For example, Stock and Watson (2002) identify predictive factors for key macroeconomic variables using linear regressions. In finance, a prominent application is the prediction of equity premia using financial and economic variables, as examined by Welch and Goyal (2008).

A critical challenge for such models in the time series setting is that their coefficients often appear to undergo structural changes due to shocks such as policy shifts, technological advances, or evolving consumer and investor behaviour. Model instability can undermine both in-sample fit and out-of-sample performance. Detecting change points in linear models is hence often a critical first step toward using them in practice. Most existing detection methods assume stationary and homoscedastic error terms; see Chapter 4 of Horv\'ath and Rice (2024), Chapter 4 of Chen and Gupta (2014), and Niu et al. (2016) for a review of change point detection methods for linear models. The assumption of homoscedasticity often appears to be implausible in practice, as model residuals frequently exhibit heteroscedasticity as well as changes in their distribution coinciding with other changes in the model parameters. This paper focuses on adapting stability tests for linear models to accommodate heteroscedastic covariates and errors.

The effect of heteroscedasticity in change point analysis has drawn increasing attention recently. Zhou (2013) and Xu (2015) advise that commonly used CUSUM-based change point procedures can become over-sized and unreliable in the presence of change points in the variance of the error process. To deal with this issue, several methods have been proposed to adapt limits for classical CUSUM-type statistics under heteroscedasticity. In the setting of changes in the mean of scalar time series, Zhou (2013) suggests a wild-bootstrap procedure to estimate the limiting distribution. Astill et al. (2023) develop a CUSUM based monitoring scheme for financial data allowing for time varying volatility. Xu (2015) builds a time transformed Wiener process, and G\'{o}recki et al. (2018) make use of Karhunen--Lo\'{e}ve expansions to characterize the limit of statistics based on heteroscedastic observations. Horv\'ath et al. (2021) derive a Wiener process-based limit for heavily-weighted CUSUM processes constructed from linear model residuals. Georgiev et al. (2018) consider the change point detection problem in predictive regression models allowing for non-stationary covariates.

In this paper, we consider a linear regression model for a scalar response $y_i$ on a $d$-dimensional covariate ${\bf x}_i$, with $R$ possible changes:

where $({\bf x}_1,y_1),...,({\bf x}_N,y_N)$ are the observed data, ${\bf x}_i \in \mathbb{R}^d$ and $y_i \in \mathbb{R}$. The regression parameter changes from $\boldsymbol \beta_{\ell}$ to $\boldsymbol \beta_{\ell+1}$ at the potential change points $k_1, \ldots, k_R$. When for example the covariates contain lagged values of an exogenous series or the response, (ref) becomes a predictive regression model with changing coefficients. We are interested in testing the null hypothesis that the regression parameter remains constant over the sample period:

versus the alternative hypothesis that there exists at least one change point,

Under $H_0$ we denote the common regression parameter as $\boldsymbol \beta_0$, which can be estimated by the least squares estimator $$ \hat{\boldsymbol \beta}_{N}=\left( {\bf X}_{N}^\top{\bf X}_{N}\right)^{-1}{\bf X}_{N}^\top{\bf Y}_{N}, $$ where ${\bf Y}_{N}=(y_1, \ldots, y_N)^\top$ is the vector containing the responses, and the design matrix is given by ${\bf X}_N = ({\bf x}_1 \mid \cdots \mid {\bf x}_N)^\top$. Thus, the linear model residuals are computed as $$ \hat{\epsilon}_i=y_i-{\bf x}_i^\top\hat{\boldsymbol \beta}_N, \quad 1\leq i \leq N. $$ The maximally selected F--tests of ${H}_0$ may be expressed as functionals of the standard CUSUM process of the covariate weighted residuals

where $\underset{\emptyset}{\sum}=0$. It follows that ${\bf Z}_N(t)=\mathbf{0}$, if $t\in [0,1/(N+1))$ and $t\in(N/(N+1),1]$. Most existing methods to test for change points in linear models make use of functionals of ${\bf Z}_N$. The likelihood ratio based tests proposed by Bai (1995, 1997a,b, 1998) and Bai and Perron (1998, 2003) are asymptotically equivalent with functionals of ${\bf Z}_N$. For example, the likelihood ratio based method in Bai and Perron (1998) can be written as the maximum of the standardized increments of the process ${\bf Z}_N$. Similarly, Bai (1999) develops a maximally selected least squares test to determine whether $R$ or $R+1$ changes are present in model (ref). It is also asymptotically equivalent with a functional of ${\bf Z}_N$. Hidalgo and Seo (2013) considers Lagrange multiplier and maximum likelihood statistics, respectively, for testing the constancy of parameters in parametric time series models, which are also asymptotically equivalent to functionals of ${\bf Z}_N$ under model (ref).

We provide in this paper a comprehensive asymptotic analysis under $H_0$ of the weighted functionals of ${\bf Z}_N$ allowing for quite general forms of heteroscedasticity in the covariates and the errors in (ref). In particular, we consider a model for non-stationary errors allowing for both {\it smooth} and {\it abrupt} changes in the error variance. An interesting consequence of the results presented is that asymptotics for the CUSUM process of the unobservable series $\{{\bf x}_i\epsilon_i\}$ and for the observable series $\{{\bf x}_i\hat{\epsilon}_i\}$ are the same with homoscedastic covariates/errors, although this does not remain true in heteroscedastic scenarios. If the volatility of the covariates ${\bf x}_i$ changes during the observation period, then the asymptotic distribution of the weighted CUSUM is affected by the estimation of the regression parameter. However, the asymptotic results for suitably standardized CUSUM statistics will, interestingly, still satisfy Darling--Erd\H{o}s type limit results in this case. The behaviour of these statistics under $H_A$ is also detailed, and it is shown that the typically weighted functionals of ${\bf Z}_N$ are consistent in detecting multiple change points.

The finite sample performances of the proposed tests are compared and studied in a Monte Carlo simulation study, which supports that the adaptations proposed to handle heteroscedasticity of the errors and to improve finite sample performance work well in practice and outperform the existing approaches of Xu (2015), Perron et al. (2020) and Horv\'ath et al. (2021). We then illustrate the proposed methods through an application to testing model instability in macroeconomic variables and equity return prediction models.

The rest of the article is organized as follows. In Section (ref), we detail the asymptotic theory for several commonly used functionals of ${\bf Z}_N$. Section (ref) extends the results for a model with more generally non-stationary errors. Section (ref) details the computation of critical values and assesses the finite-sample performance of the proposed tests through Monte Carlo simulations, comparing them with existing methods. Data applications are given in Section (ref), and Section (ref) concludes with some remarks.

Abrupt changes in the variance model

Let $\{{\bf z}_i=({\bf x}_i^\top, \epsilon_i)^\top, -\infty<i<\infty\}$ denote the process describing the covariates and error terms in (ref). We first consider the case in which the second order properties of ${\bf z}_i=({\bf x}_i, \epsilon_i)^\top, 1\leq i \leq N$ may change $M$ times during the observation period, where $1 < m_1<m_2<\ldots<m_M < N$ denote the times at which the second order properties of the ${\bf z}_i$'s might change. We assume

Leybourne et al.\ (2006), Pein et al.\ (2017) and Horv\'ath et al.\ (2021) introduce heteroscedastic models, similar to Assumption (ref), in change point analysis. We use a decomposable Bernoulli shift model for each segment of stationarity. Let $m_0=0, m_{M+1}=N$ and correspondingly $\tau_0=0$ and $\tau_{M+1}=1$. We note that the unknown times $m_i$ may or may not coincide with the change point locations $k_j$.

{Below $\left\| \cdot \right\|$ denotes the Euclidean norm.}

Under Assumption (ref), the errors and covariates are not stationary over the whole observation period, but are drawn from a stationary process on the sub-segments $(m_{\ell-1}, m_\ell]$, $1 \leq \ell \leq M+1$. Conventionally, the error terms $\epsilon_i$ are independent of ${\bf x}_i$, and are homoscedastic in the sense that the variance of the conditional distribution of $\epsilon_i$ given ${\bf x}_i$ remains constant with respect to $i$. This condition might not hold under Assumption (ref), resulting in a heteroscedastic model. In Section (ref), we further relax this assumption to allow for non-stationarity within each sub-segment.

In this section, we aim to establish the asymptotic behaviour of ${\bf Z}_N(t)$, as defined in (ref), under the null hypothesis of no change in the regression parameter under Assumption (ref).

In order to identify the regression parameters, we require

Assumption (ref) postulates that the identification of the regression parameters holds on all sub-segments of stationarity. To state the weak limit of the process ${\bf Z}_N$, we also need to introduce $M+1$ long run covariance matrices reflecting the changing covariances between stationary subintervals. Let

We show that these matrices are well defined under Assumption (ref). We define the process $\{\boldsymbol \Gamma(t), 0\leq t \leq 1\}$ as

where $\{ {\bf W}_{{\bf D}_j}(t), \; t \ge 0\}, 1\leq j\leq M+1$, are independent $d$ dimensional Brownian motions such that $E{\bf W}_{{\bf D}_j}(t)=\bf0$ and $E{\bf W}_{{\bf D}_j}(t){\bf W}_{{\bf D}_j}^\top(s)=\min(t,s){\bf D}_j$. The process $\boldsymbol \Gamma(t)\in \mathbb{R}^d$ is Gaussian process with $E\boldsymbol \Gamma(t)=\bf0$,

and

As it is common in the theory of linear regression, we require

Note that under Assumption (ref), the matrices ${\bf A}_i$ are well defined. Let

where

for $\tau_{k-1}< t \leq \tau_k, 1\leq k \leq M+1$.

One often considers weighted functionals of ${\bf Z}_N$ to improve the power of tests, particularly when changes occur near the ends of the sample. In our context, we apply the weight function $w(t)$ satisfying the following properties:

Due to using the weight function $w$, we need an integral condition for the existence of a non-degenerate limit distribution of the weighted process ${\bf Z}_N$. Let

The integral $I(w,c)$ characterizes the upper and lower classes for the Brownian bridge at 0 and 1 (It\^o and McKean, 1965; O'Reilly, 1974). According to It\^{o} and McKean (1965), we know that $I(w, c)< \infty$ with some $c$ if and only if

where $\{B(t), 0\leq t \leq 1\} $ is a Brownian bridge. The most commonly used weight function satisfying $I(w,c)<\infty$ and Assumption (ref) is $w(t)=[t(1-t)]^\kappa$, $0\leq \kappa <1/2$. Since $I([t(1-t)]^{1/2},c)=\infty$ for all $c>0$, (ref) cannot hold with $w(t)=[t(1-t)]^{1/2}$, so in this case we have a different limit distribution. We note that applying the weight function $w(t) = [t(1 - t)]^{1/2}$ can lead to a slow convergence rate. We refer to Cs\"org\H{o} and Horv\'ath (1993) for more details and discussions on weighted empirical and Gaussian processes.

We thus have $E\bar{\boldsymbol \Gamma}(t)=\bf0$ and

We note if $\{{\bf x}_i, -\infty <i<\infty \}$ is stationary, then ${\bf A}_1=\dots={\bf A}_{M+1}$ and therefore ${\bf u}(t)=\mathbf{O}$. In this case

{ It is important to note that due to the estimation of $\boldsymbol \beta$, the weighted CUSUM's of ${\bf x}_i\epsilon_i$ and ${\bf x}_i\hat{\epsilon}_i$, $1\leq i \leq N$, have a qualitatively different asymptotic distribution in heteroscedastic models. If for instance ${\bf A}_1 = \dots = {\bf A}_{M+1}$, i.e., the sequence $\{{\bf x}_i, -\infty <i<\infty\}$ is homoscedastic, and the heteroscedasticity is only in the errors $\{\epsilon_i, 1\leq i \leq N\}$, then ${\bf v}(t)=\mathbf{O}$, so the effect of estimating $\boldsymbol \beta$ does not appear in the limit distribution. In other words, interestingly the limit behaviour of the CUSUM of ${\bf x}_i\epsilon_i$'s and ${\bf x}_i\hat{\epsilon}_i$'s are the same if $\{{\bf x}_i, -\infty <i<\infty\}$ is homoscedastic even if $\{\epsilon_i, -\infty <i<\infty\}$ is not.}

We now turn our focus to the standardized statistics using the weight function $w(t)=[t(1-t)]^{1/2}$, and we use the following assumption

Assumption (ref) is a sufficient condition to obtain limit results for suitable standardized weighted supremum functionals of ${\bf Z}_N$. The proofs in Appendix (ref) show that under $H_0$ the maximum of such functionals is asymptotically reached on the intervals $(0, \tau_1]$ or $(\tau_M, 1]$, and therefore the maximum taken on these intervals determines the limit distribution. We can standardize the CUSUM process with the matrix valued function $$ \tilde{{\bf G}}(t)=E\bar{\boldsymbol \Gamma}(t)\bar{\boldsymbol \Gamma}^\top(t),\quad 0<t<1. $$ Let us consider the following “Darling--Erd\H{o}s" type statistics, which converge weakly to extreme value laws.

We note the weight function $1/w(t)$ is not explicitly used in the above statistics, as the statistics have been inherently weighted through the standardization term $\tilde{{\bf G}}(t)$.

The matrix valued function $\tilde{{\bf G}}(t)$ is unknown, and in practice must be estimated from the sample. In order to do so, first we estimate ${\bf G}(u)$ with a long run variance kernel estimator based on a fraction $u$ of the data, where $0<u\leq 1$. Here we use the standard kernel covariance estimator for the weighted residuals $\{ {\bf x}_i\hat{\epsilon}_i, 1\leq i \leq N \}$. Considering the kernel function $K$ and a bandwidth parameter $h=h(N)$, we require

The parameter $\rho$ in Assumption (ref)(v) indicates the order of the kernel near zero, which also approximates the asymptotic bias of kernel--based long run variance estimators. For example, the popular Bartlett kernel has order $\rho=1$ and the Parzen kernel has order $\rho= 2$ (see e.g. Andrews, 1991). We then have the long run covariance matrix estimator ${\bf G}_N(u)$ that is computed from the weighted residuals $\{ {\bf x}_i\hat{\epsilon}_i, 1\leq i \leq \lfloor Nu\rfloor\}$. Let

Now $\{{\bf G}_N(u), 0\leq u \leq 1\}$ is defined as

where $K$ and $h$ satisfy Assumption (ref).

To approximate the limit in Theorem (ref), and according to (ref), we need to estimate ${\bf u}(t)$. We use

Then, we can estimate $\bar{{\bf G}}(t,s)$ to approximate the limit in Theorem (ref) using the plug-in estimator by replacing ${\bf G}(t)$ with ${\bf G}_N(t)$ and ${\bf u}(t)$ with ${\bf u}_N(t)$ in (ref). For $\tilde{{\bf G}}(t)$ in Theorem (ref), we can use the plug--in estimator

The consistency of $\tilde{{\bf G}}_N(t)$ follows from Theorem (ref).

We now turn to establishing the behavior of the functionals of ${\bf Z}_n$ under the alternative hypothesis. To do so, we further assume that

Under Assumption (ref), the potential change points are well separated. Also, changes in the parameters can degenerate along with the sample size increasing, i.e.,\ it is possible to have $\|\boldsymbol \beta_\ell-\boldsymbol \beta_{\ell-1}\|=o(1)$.

Let

where $| \cdot | $ is the Lebesgue measure. Next we define

if $\theta_{k-1}<t\leq \theta_k, 1\leq k \leq R+1$.

If the changes in the volatility occur at the same time as in the linear model coefficients, i.e., $M=R$, and $r_\ell=m_\ell$ for all $1\leq \ell \leq R$, then the formulas for $\boldsymbol \beta^{**}$ on $\bar{{\bf g}}_N(t)$ are simpler. In this case $$ \boldsymbol \beta^{**}=\left( \sum_{j=1}^{R+1} (\theta_j-\theta_{j-1}) {\bf A}_j \right)^{-1} \sum_{\ell=1}^{R+1}(\theta_\ell - \theta_{\ell-1}) {\bf A}_\ell \boldsymbol \beta_\ell, $$ and $$ \bar{{\bf g}}_N(t) = \sum_{j=1}^{k-1} (\theta_j - \theta_{j-1}){\bf A}_j(\boldsymbol \beta_j-\boldsymbol \beta_j^{**})- t\sum_{\ell=1}^{R+1}(\theta_\ell-\theta_{\ell-1}){\bf A}_\ell(\boldsymbol \beta_\ell - \boldsymbol \beta^{**}), $$ for $\theta_{k-1}<t\leq \theta_k$, $1\leq k \leq R+1$.

The tests will stay consistent if $\tilde{{\bf G}}(t)$ is replaced $\tilde{{\bf G}}_N(t)$ in the testing procedures. For example, if $\underset{N\rightarrow \infty}{\lim}(N/h)^{1/2}(\log \log N)^{-1/2}\underset{1\leq \ell \leq M+1}{\max} \| \boldsymbol \beta^{**}-\boldsymbol \beta_\ell\|\rightarrow 0$, then $$ \left(\log \log N \right)^{-1/2} \underset{0<t<1}{\sup} \frac{1}{(t(1-t))^{1/2}} \left( {\bf Z}_N^\top(t) \tilde{{\bf G}}_N^{-1}(t) {\bf Z}_N(t) \right) \overset{P}{\rightarrow} \infty. $$ The proofs of the theorems in this section are provided in the Online Supplement Section (ref).

Smoothly changing error variance model

So far, based on Assumption (ref), we have considered a model where the structure of the errors in the observations might change during the observation period, but that the errors are piecewise stationary. In this section, we extend to the case when there are smooth changes in the variance of the errors during the intervals $(m_{i-1}, m_i], 1\leq i\leq M+1$. Inspired by the mean change point model in G\'orecki et al.\ (2018), we modify the model of (ref) as

where the errors $\{\epsilon_i, -\infty<i<\infty\}$ are as described in Section (ref), and it holds that

Hence, we allow the variances of the errors $g(i/N)\epsilon_i$ to change even within the intervals $m_{\ell-1} < i \leq m_\ell$, for $1 \leq \ell \leq M+1$. The asymptotic behavior of ${\bf Z}_N$, for $0 \leq t \leq 1$, in model (ref) is established below.

We note that $\mbox{\boldmath${\Lambda}$}(t)$ is a $d$ dimenional time transformed Brownian motion with $E\mbox{\boldmath${\Lambda}$}(t)=0$, and $$ E\mbox{\boldmath${\Lambda}$}(t)\mbox{\boldmath${\Lambda}$}^\top(s)={\bf {H}}(\min(t,s)), $$ with

Let $$ \mathcal {a}(t)=\int_0^tg^2(u)du. $$ We note that the variance of the coordinates of $\mbox{\boldmath${\Lambda}$}(t)$ are proportional to $\mathcal {a}(t)$, so it is natural to assume

We now present the weighted version of Theorem (ref).

Theorem (ref) implies immediately that

and

The proof of Theorem (ref) is given in the Supplemental Section (ref). Since we modify only the error term, Theorem (ref) also holds under model (ref).

Next, we consider the standardized statistics of Theorem (ref) when

Let $$ \bar{{\bf {H}}}(t)={\bf {H}}(Nt)-2t{\bf {H}}(tN)+t^2{\bf {H}}(1). $$ We note that $\bar{{\bf {H}}}(k/N)$ is the covariance matrix of $\mbox{\boldmath${\Lambda}$}(k)-(k/N)\mbox{\boldmath${\Lambda}$}(N)$.

The main difference between Theorems (ref) and (ref) arises from the substantially smaller variances of $g(i/N)\varepsilon_i$ under (ref) when $i$ is small, compared to the variances when $i$ is close to $N$. The variances of $g(i/N)\epsilon_i$ are converging to 0 as $i/N\to 0$, while to $g^2(1) E\epsilon^2_0$, if $i/N\to 1$. We then show that Theorem (ref) remains true if the variances of $g(i/N)\epsilon_{i}$ converge to positive constants if $i/N\to 0$ and $i/N\to 1$. We replace (ref) with

We also note that

under the assumptions of Theorem (ref). This means that there is no difference between the covariance matrix estimated Darling--Erd\H{o}s results in Sections (ref) and (ref). No information on the error structure is required to implement the testing procedure, as long as $E\varepsilon_i^2 \geq c > 0$. We can use the same method to compute the critical values used in cases of abrupt change errors or smoothly changing variance errors.

Monte Carlo simulation and Finite Sample Performance

The computation of critical values

\\ To assess the finite-sample performance of our tests, we focus on the standardized Darling--Erd\H{o}s-type statistics presented in Theorems (ref) or (ref). As noted in Remark (ref), the standardized statistics are more practical in data applications, as they do not require prior information on heteroscedasticity. However, Darling--Erd\H{o}s-type statistics typically suffer from slow convergence due to their exponential-type limiting distribution, often leading to tests that are under-sized and with reduced power. We put forward an improved finite sample approximation of the Darling--Erd\H{o}s limiting distribution in this section. Let

where $\tilde{{\bf G}}_N(t), 0<t<1$ is defined in (ref). The long run variance estimators ${\bf G}_N(u)$ in (ref) were computed with the Bartlett kernel, and the bandwidth is selected through the automatic bandwidth selection method of Andrews (1991). This subsection explains the computation of critical values for the test statistics $V^{HET}_N(1/2)$ and $Q^{HET}_N(1/2)$.

In Appendix (ref), we show that in Theorem (ref), a Gaussian approximation is first established for the process ${\bf Z}_N$, and then the limiting distribution for the maximum of Gaussian processes is applied. The Gaussian approximation yields specifically that

for any $c_1(N),c_2(N) \rightarrow 0$ and satisfying $c_1(N)= O\left((\log N)^{\kappa_1}/N\right)$ and $c_2(N)= O\left( (\log N)^{\kappa_2}/N \right)$ with any $-\infty < \kappa_1, \kappa_2< \infty$. Similarly,

where $\{B_1(t), 0\leq t \leq 1\}, \ldots, \{B_d(t), 0\leq t \leq 1\}$ are independent Brownian bridges.

To obtain the critical values, in practice one can simulate the random variable, $$ \sup_{c_1(N)\leq t \leq 1-c_2(N)} \frac{1}{(t(1-t))^{1/2}}\left( \sum_{i=1}^{d} B_i^2(t)\right)^{1/2} $$ for choices $c_1(N)$ and $c_2(N)$ such as $c_1(N)=c_2(N)=1/N$, and take its quantiles as critical values. We now instead introduce an approximation inspired by Vostrikova (1981), which does not require any simulation to obtain the empirical quantiles of the limits. We recall from Cs\"org\H{o} and Horv\'ath (1997, p.\ 366),

where $$ U^*_d(x)=\left(\sum_{i=1}^d U_i^2(x)\right)^{1/2}, $$ and $\{U_1(x), -\infty <x <\infty\}, \ldots, \{U_d(x), -\infty<x<\infty\}$ are independent, identically distributed Ornstein--Uhlenbeck processes, i.e.\ Gaussian processes with $EU_i(x)=0$ and $EU_i(x)U_i(y)=\exp(-|x-y|/2)$. As a result, we get from (ref) and (ref) that

and

The critical values of the statistics $V^{HET}_N(1/2)$, and $Q^{HET}_N(1/2)$ are be approximated by using (ref) and (ref). In Vostrikova (1981), it is shown that for $T>0$ and $r\geq 1$,

where $\Gamma(\cdot)$ denotes the Gamma function. We ignore the $O(1/x^4)$ term, and then compute critical values directly from (ref).

In an unreported simulation, we found that using critical values based on the Vostrikova approximation yields higher power than those from the Darling--Erd\H{o}s limit. Therefore, we use the Vostrikova-based critical values in the analysis below. For further details on the implementation of the non-standardised CUSUM statistics in Theorems (ref), (ref), and (ref), we refer the reader to Section (ref) of the Supplementary Material.

Monte Carlo Simulations

\\ We now introduce the model and settings for the Monte Carlo simulation study, with results discussed in the next subsection. We consider a data generating process (DGP) taking the form (ref) and simply allow one change point in the regression parameters:

The covariates are ${{\bf x}}_i = (1, x_{2,i})^\top$ for $1 \leq i \leq N$, where $x_{2,i}$ is a heteroscedastic covariate following a segmented independent and identically distributed (i.i.d.) normal distribution with mean 0, standard deviation 3 before, and 0.5 after the break at $\lfloor 0.45 N \rfloor$. Prior to the change point $k_1$, we set the regression parameter as $\boldsymbol{\beta}_1 = (1 , 1 )^\top$, and the coefficients become to $\boldsymbol{\beta}_2 = (1+\delta, 1+\delta)^\top$ after the change. The change point size $\delta$ varies in the range $\delta \in \{-1.5, -1.2, -0.9, -0.6, -0.3, 0, 0.3, 0.6, 0.9, 1.2, 1.5\}$. The null hypothesis $H_0$ holds when $\delta=0$. The rejection rates at the nominal size $\omega = 0.05$ are reported as power curves in terms of the change size $\delta$, for each DGP and for middle change point at $k^ = \lfloor 0.5N \rfloor$ and early change point $k^ = \lfloor 0.2N \rfloor$.

For the error term $\epsilon_i$, we consider four heteroscedastic processes, encompassing both abrupt variance changes (cases i--iii) and a smooth variance change (case iv).

(i) (Normal) the error terms $\epsilon_i$ are i.i.d normal random variables:

(ii) (AR) the error terms $\epsilon_i$ follow autoregressive (AR-1) process:

(iii) (GARCH) the error terms $\epsilon_i$ follow a stationary GARCH(1,1) process defined by

and the $\varepsilon_i$'s are i.i.d. standard normal random variables. In this exercise, we always set $m^*=\lfloor 0.45N \rfloor$, i.e., there is a variance change around the middle of the sample.

(iv) (HeteSmooth) Lastly, the error term follows a heteroscedastic smooth variance change process,

We set the sample size $N= \{125,250\}$. All reported results are based on 2000 independent replications in each setting.

We compare the proposed methods to three existing tests in the literature that considered change point detection in linear models with heteroscedastic errors, including Xu (2015), Perron et al. (2020) and Horv\'ath et al., (2021). Xu (2015) propose CUSUM--type tests to detect changes in linear models with nonstationary variance. We adopt the robust CUSUM test, with critical values derived from the limit approximated by a T-discrete steps Wiener process, denoted as the $\rm XUC$ test. Perron et al. (2020) introduce a likelihood ratio--type test to detect changes in the coefficients, accommodating scenarios where the coefficient change up to $p$ times and the variance of the errors changes up to $q$ times. Specifically, we adopt their simulation setup, setting $p\in[0, 1, 2]$ and $q\in [1,2,3]$. The null hypothesis of no change is rejected when the test statistics exceed the critical values in at least one model specification. We employ the statistic $LR_{3,N}$ and refer to it as the $\rm{PYZ}$ test. Horv\'ath et al. (2021) introduce a R\'enyi--type statistic to compare the least squares estimators from sub-segments split by potential change points. Following their suggestion, we choose the tuning parameters $a_N = b_N = N^{1/2}$ and denote the statistic as $\rm{HMR}$ test. In Online Supplement Section (ref), we examine the effect of the choice of $\kappa$ and the methods used to calculate the critical values. The results generally recommend using Darling--Erd\H{o}s type statistics with the distributional approximation obtained from (ref) in terms of testing power. Therefore, we use the tests $V_N^{HET}(1/2)$ and $Q_N^{HET}(1/2)$ in this experiment.

Figure (ref) and (ref) display the power curves of the test candidates for $N = 125$ and $N = 250$, respectively. The overall patterns are consistent across both early and mid-sample changes, with greater power observed as the sample size increases. Tests $\mathrm{XUC}$ and $V^{\mathrm{HET}}_N(1/2)$ maintain approximately nominal size, while $Q^{\mathrm{HET}}_N(1/2)$ is slightly oversized. The $\mathrm{HMR}$ test performs well under GARCH errors but exhibits substantial size distortion when applied to other error structures, particularly under HeteSmooth errors. The test $\rm{PYZ}$ is oversized in all DGPs, and it deteriorates in models generated with GARCH errors. We note that the test $\rm{PYZ}$ is designed to mitigate power reduction, but the simulation shows that it can be oversized when dealing with changes on the tails and heteroscedastic covariates that fall outside the scope of its intended design.

Under the alternative hypothesis, all tests start to gain reasonable power in large samples. The test $\rm{XUC}$ exhibits the lowest overall power, particularly with small sample sizes and early change points. The test $\rm{HMR}$ is relatively competitive in models with NORMAL and GARCH errors. The test $\rm{PYZ}$ achieves significantly enhanced power, but it is not reliable due to the size distortions. The test $Q^{HET}_N(1/2)$ consistently outperforms $V^{HET}_N(1/2)$, which is primarily used in Section (ref) for empirical applications. Overall, the simulation results imply distortions of oversize or reduced power of the existing tests by encountering severe heteroscedasticity and various locations of changes in the linear models.

Additional simulation results for models with homoscedastic errors, as well as the effects of change-point locations and signal-to-noise ratios, are reported in Online Supplement Section (ref).

Empirical data examples

Testing instability in macroeconomic prediction

\\ We first illustrate a data application of the proposed methods to test the instability of predictive regression models. A typical univariate predictive regression model takes the form

where $x_i$ denotes the predictor variable, typically representing observation from lagged values. These simple models are widely used in macroeconomic studies, which we also analyze here.

We follow McCracken and Ng (2016) and aim to forecast monthly U.S. GDP growth, industrial production, nonfarm employment, and total CPI inflation, indexed by $j = \{1, 2, 3, 4\}$. The covariate that we use to forecast each of these series is the real activity/employment predictive factor derived from a panel of 134 U.S. macroeconomic indicators in McCracken and Ng (2016). We consider a univariate predictive model because, as stated in McCracken and Ng (2016), two factor models only provide slight improvements marginally in terms of forecasting error. The sample ranges from January 1993 to December 2022 with 357 observations in total\footnote{The monthly macroeconomic data is collected from the webpage “https://research.stlouisfed.org/econ\\/mccracken/fred-databases/".}. {Preliminary data analysis showed that the covariates appeared to exhibit abrupt heteroscedasticity at the onset of several periods of market turbulence, such as the Great Recession of 2008-2009 and the beginning of the COVID-19 pandemic. The standard deviation of the model residuals also shows changes when using a rolling window estimation, see Figure (ref).}

The proposed tests were applied to each model for detecting change points. A change point is estimated based on each test statistic using the argument at which the corresponding normalized CUSUM processes achieved their maxima. We apply standard binary segmentation based on each change point test statistic with a threshold taken to be the 95% null significance level in an attempt to detect additional change points. Table (ref) in the appendix shows the change points detected by each approach and the corresponding estimated model coefficients. Consistent with the simulation results, the statistics $Q^{HET}_N(1/2)$ detects more changes, while the test $V^{HET}_N(1/2)$ is relatively conservative.

Table (ref) shows the coefficient estimations in subsamples split by the first four estimated change points. In general, we found changes occurring in the GDP growth, industrial production and nonfarm employment models. Indicated by the $V^{HET}_N(1/2)$ and $Q^{HET}_N(1/2)$ tests, these three models experience changes around the period of recovery from the 2008 great recession. The tests also suggest changes around the COVID-19 pandemic in GDP growth and nonfarm employment models, with $Q^{HET}_N(1/2)$ detecting additional changes during the 1997 Asian financial crisis in the GDP growth model and during the early 2000s recession in all three models.

Figure (ref) illustrates an example of segmentation produced by $Q_N^{HET}(1/2)$ method for predictive regression model for nonfarm employment in terms real activity/employment predictive factor. The black line shows the model residuals from the model (ref) applied to the entire sample, and the blue line shows rolling window estimates of the standard deviation of the model residuals, with window size 24 months. The shaded bars indicate the U.S. business cycle contractions according to NBER\footnote{NBER U.S. business cycle dating: “https://www.nber.org/research/business-cycle-dating".}.

Changes in investor sentiment effect on the U.S. stock market

\\ In this section, we demonstrate a second application to detect changes in explaining the sentiment anomaly on cross--sectional U.S. stock returns. This provides an example of the proposed tests for an explanatory regression model with multiple covariates. The aims of modeling market sentiment anomalies are to try and model two phenomena; how the demand for speculative investments drives stock prices away from their fundamental values, and relative arbitrage, i.e.\ the existence of a collection of stocks that are too risky and costly for arbitrage. The topic-influential work of Baker and Wurgler (2006) reviews the anecdotal history of investment sentiment in the U.S.\ between 1961 and 2002, and constructs sentiment factors to predict cross--sectional stock returns.

We procured a dataset covering the period January 1960--June 2022\footnote{The data is collected from the webpage of Jeffrey Wurgler “http://people.stern.nyu.edu/jwurgler/".}, containing two sentiment factors (SENTI) constructed from six underlying sentiment proxies, including the closed--end fund discount, NYSE share turnover, the number and average first--day returns on IPOs, the equity share in new issues, and the dividend premium. Our study uses the second sentiment factor because it accounts for business cycle variations. The sample consists of 684 time series observations, and extends beyond the data considered in Baker and Wurgler (2006); it includes some recent major economic events of note, such as the US housing bubble, the great recession, and the Covid-19 pandemic.

Following Baker and Wurgler (2006), we consider the premium of the size factor (small--minus--big, SMB) factor as a dependent variable to verify the distinct sentiment effect among small firms. The linear regression model specifies four dependent variables, including the sentiment index, market excess return (RMRF), premium of the book--to--market (high--minus-low, HML), and premium on winners minus losers (momentum, MOM) factors\footnote{The Fama--French--Carhart factors are obtained from the data library “https://mba.tuck.dartmouth.edu\\/pages/faculty/ken.french/data_library.html".}. The model can be stated explicitly as:

{ While all other covariates appear to be reasonably stationary, the covariate SENTI exhibited apparent changes in its variance. The standard deviation of the model residuals again exhibits heteroscedasticity when using a rolling window estimation. }

We test for change points in the regression parameters in (ref) using the statistics $V^{HET}_N(1/2)$ and $Q^{HET}_N(1/2)$. The changes around April 2002, April 2009, and June 2020 are detected for both tests in Table (ref), indicating changes occurring during the burst of the early 2000s recession, the recovery from the great recession, and the outbreak of the Covid-19 pandemic. The $Q^{HET}_N(1/2)$ test suggests the presence of one more change in July 1973, which can be a consequence of the 1973–-1974 stock market crash.

Table (ref) also shows the results of the estimations when the sample is segmented with the detected changes. Focusing on $V^{HET}_N(1/2)$ and $Q^{HET}_N(1/2)$, we estimate $\beta_2=-0.22$ with p--value 0.07 based on the sampling period between 1965 and 2002, while $\beta_2$ is estimated as $-0.23$ and $-0.18$ in subsamples respectively by split from July 1973. Both coefficients appear to be insignificant, but their p--values close to 0.10. Our findings are roughly consistent with Baker and Wurgler (2006), who estimated $\beta_2=-0.30$ with p--value $0.15$ using the period 1961 to 2002. The negative sign of the coefficient indicates that there is a negative relationship with the sentiment premium, i.e., the small firms turn to gain less returns with intense market sentiment. This effect becomes more manifest during the formation and collapse of the U.S. housing bubble between 2002 and 2009, given the coefficient $\beta_2$ enlarges to $-1.29$. The sentiment effect then becomes insignificant in subsamples after 2009, but the uncovered negative effect is consistent throughout each subsamples.

Conclusion

We propose quadratic forms and maxima of weighted CUSUM residual processes to test for multiple changes in linear model parameters under potential heteroscedasticity in both covariates and errors. The error variance is allowed to change either abruptly or smoothly. The asymptotic distributions of the proposed test statistics are established under general conditions that accommodate both homoscedastic and heteroscedastic cases. We examine the finite sample performance of the standardized statistics in detail. Monte Carlo simulations demonstrate that the tests exhibit good size and power in finite samples, and that the adjustments for heteroscedasticity in model errors perform well in practice. We applied our method to find changes in popular macroeconomic and return prediction models, and to detect changes in the sentiment asset pricing models in the U.S.\ stock market.

\setcounter{section}{0}