Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
87,064 characters · 9 sections · 98 citation commands
Detecting Multiple Structural Breaks in Systems of Linear Regression Equations with Integrated and Stationary Regressors
\sloppy
\singlespacing
\thispagestyle{empty}
\onehalfspacing
Accounting for structural breaks is crucial in time series analysis, particularly in settings involving long spans of data, where the models are more likely to be affected by multiple structural breaks. More specifically, we focus on systems of equations with a mix of integrated and stationary regressors. Thus far, the literature on structural breaks has provided only few methods applicable to linear regressions with multiple equations and integrated regressors BaiLumsdaineStock1998, LiPerron2017, OkaPerron2018. Without prior knowledge about the structural breaks, methods are needed that precisely determine the number of structural breaks, their timing, and simultaneously estimate the model's coefficients.
Considering the increasingly larger sample sizes in empirical studies, MacKinnon2023 cautions against the use of algorithms that are $O(T^{\eta})$ for $\eta \gg 1$ which become impractical for sufficiently large $T$. The currently available methods to solve change-point problems in models that we consider in this paper consist of likelihood-based approaches using dynamic programming techniques BaiPerron1998, BaiPerron2003 and are characterized by computational costs quadratic in the number of observations. While the approaches are generally very precise, they are not computationally efficient in situations where $T$ is large. Consequently, we propose to lower the computational burden and design a a structural break detection algorithm that is feasible for large $T$ systems.
We consider a penalized regression approach based on the group LASSO estimator to account for multiple structural breaks in such systems which, to the best of our knowledge, has not been explored in the literature yet. Although estimators based on the penalized regression principle have become popular in the context of change-point problems, few prior studies apply them to linear regressions with integrated regressors SchmidtSchweikert2019, Schweikert2021 or linear regressions with multivariate responses Gao2019, SafikhaniShojaie2020. While existing approaches follow a specific-to-general principle utilizing a likelihood-based approach to sequentially increase the number of breaks in a model BaiLumsdaineStock1998, QuPerron2007, LiPerron2017, OkaPerron2018, we take a general-to-specific approach shrinking down the number of breakpoint candidates to find the best fitting model. While the likelihood-based approach employs dynamic programming techniques and is computationally efficient in rather short samples with (possibly) many structural breaks, the proposed model selection approach is particularly useful for long samples with a moderate number of structural breaks. It can be shown that it has computational costs linear in the number of observations. Therefore, it is a well-suited solution to account for structural breaks in the long-run relationships between trending variables.
In this paper, we extend the two-step estimator proposed in ChanYauZhang2014 for univariate structural break autoregressive (SBAR) processes. To do so, we modify the group LARS algorithm specifically tailored for univariate change-point problems and extend it to cover multivariate systems. Moreover, we generalize the model specification and allow for a mix of stationary and integrated regressors as well as deterministic trends. Consequently, our approach is flexible enough to model structural breaks in several special cases like, for example, seemingly unrelated regression (SUR) models and dynamically augmented cointegrating regressions.\footnote{While the model structure, in principle, includes the possibility to consider piece-wise stationary VAR models, our technical analysis relies on Assumption (ref) stated below which is not compatible with VAR models. We refer to SafikhaniShojaie2020 who use a slightly different penalty to cover a high-dimensional version of this case.}
The idea to perceive the change-point problem in linear regressions as a model selection problem has spawned a diverse literature Harchaoui2010, BleakleyVert2011, ChanYauZhang2014, SafikhaniShojaie2020, Schweikert2021. In principle, it is possible to shift and turn the regression hyperplane at every point in time using appropriate indicator variables. Finding the true structural breaks corresponds to selecting relevant indicators and eliminating irrelevant indicators thereby optimizing the fit under sparsity. This leads to a high-dimensional penalized regression model with the total number of parameters of the model close to the number of observations. LASSO-type estimators, introduced by Tibshirani1996, have attractive properties in high-dimensional settings with a sparse model structure. Their objective function includes a penalty for nonzero parameters and a tuning parameter controls the sparsity of the selected model. However, quite restrictive regularity conditions about the design matrix (restricted eigenvalue condition BickelRitovTsybakov2009 or strong irrepresentable condition ZhaoYu2006) need to be imposed to ensure simultaneous variable selection and parameter estimation consistency. Unfortunately, these conditions are usually violated in change-point settings, where adjacent columns of the design matrix differ only by one entry and the design matrix is highly collinear if the sample size grows large. Consequently, the conventional LASSO-type estimators need to be improved to both estimate and select the true model consistently.
Harchaoui2010 are among the first to use penalized regression methods to detect structural breaks. They focus on a piecewise constant white noise process and detect structural breaks using a total variation penalty. BleakleyVert2011 use the group fused LASSO for detection of piecewise constant signals and ChanYauZhang2014 develop the aforementioned two-step method. Ciuperca2014, JinWuShi2016 and QianSu2016b consider LASSO-type estimators for the detection of multiple structural breaks in linear regressions. BehrendtSchweikert2020 propose an alternative strategy to eliminate superfluous breakpoints identified by the group LASSO estimator. They suggest a second step adaptive group LASSO which performs comparably to the backward elimination algorithm suggested in ChanYauZhang2014. SchmidtSchweikert2019 consider cointegration tests in the presence of structural breaks in the long-run relationship and estimate those breaks with an adaptive LASSO estimator. Schweikert2021 uses the adaptive group LASSO estimator to estimate structural breaks in single-equation cointegrating regressions. The estimator developed in this paper can be understood as an extension of the Schweikert2021 approach to multiple equation cointegrated systems, so that structural change of more than one equilibrium relationship can be modelled.
Work on multiple structural change models in the context of a system of multivariate equations is relatively scarce. Quintos1995, Quintos1997 considers a general time-varying structure for the reduced-rank matrix of a vector error correction model (VECM) so that both the cointegrating vector and the adjustment dynamics may change over time. Similarly, Seo1998 develops a test for changing cointegrating vectors and adjustment coefficients at a single unknown breakpoint. BaiLumsdaineStock1998 concentrate on dating and estimating a single structural break in vector autoregressions (VARs) and multiple equation cointegrating regressions. QuPerron2007 consider the restricted quasi-maximum likelihood estimation of and inference for multiple structural changes in a system of equations. A sequential break test can be used to determine the number of structural breaks. In related studies, EoMorley2015, LiPerron2017, and OkaPerron2018 extend the likelihood-based approach in several directions.\footnote{EoMorley2015 propose confidence sets for the timing of structural break estimation in multiple equation regression models. LiPerron2017 introduce the concept of locally ordered breaks. They model structural breaks in systems of equations with a combination of integrated and stationary regressors, dealing with situations where the breaks cannot be separated by a positive fraction of the sample size. OkaPerron2018 highlight that the estimation of common breaks allows for a more precise detection of break dates in multivariate systems. They develop common break tests for this assumption. A common break is defined as a point in time at which at least one coefficient from each equation is not restricted to be the same across two adjacent segments.}
Recently, the model selection approach has been applied by Gao2019 and SafikhaniShojaie2020 to estimate change-points in a piece-wise stationary VARs. While the former study estimates the change-points for each equation separately, thereby decomposing the problem into smaller single-equation problems, the latter uses a fused LASSO penalty to deal with high-dimensional VAR systems.
In the following, we provide a rigorous analysis of the statistical properties of the proposed two-step estimator and extensive simulation experiments to analyze its finite sample properties. We conduct our technical analysis under relatively mild assumptions about the error term process. Prior studies employing LASSO-type estimators to detect structural breaks assume (Gaussian) white noise error terms ChanYauZhang2014, Gao2019, SafikhaniShojaie2020, which can be useful to model (V)AR processes. However, this assumption is too restrictive in (multiple equations) linear regressions with integrated regressors often having serially correlated errors. Naturally, it becomes more difficult to detect structural breaks if the error term process is serially correlated. Under those assumptions, we show that our estimator is able to consistently estimate the number of structural breaks, their timing, and jointly estimates the model's coefficients.
We use simulation experiments to evaluate our new approach against existing approaches like the likelihood-based approach by, inter alia, QuPerron2007, LiPerron2017, and OkaPerron2018. It is shown that the two-step estimator has competitive finite sample properties with a slight reduction in precision, but substantially improved computational efficiency. Reducing the computational burden over the likelihood-based approach is an important advantage when large sample sizes are available and a moderate number of structural breaks is expected as is often the case in empirical applications involving trending regressors. Another advantage is the joint estimation of the number of breaks, their timing, and the model's coefficients. In the likelihood-based framework, the number of breaks has to be determined based on the evaluation of two tests with the usual implications regarding size and power.\footnote{First, a double maximum test is conducted to test whether at least one break is present, then the exact number of breaks is determined testing the hypothesis of $l$ breaks versus the alternative of $l + 1$ breaks. Naturally, the sequential test procedure requires the specification of a nominal significance level $\alpha$ which implies that also for large samples, the number of breaks is overestimated in roughly $\alpha \cdot 100$% of all cases.} In contrast, the approach taken in this paper does not rely on statistical testing, instead we determine the number of breaks as the number of nonzero groups estimated by the group LASSO estimator which is then further reduced by a second step backward elimination algorithm. While we also need to ensure that each identified regime has a sufficient number of observations to estimate the coefficient changes, we need a much smaller trimming parameter than commonly applied in the literature QuPerron2007, OkaPerron2018. Hence, the two-step estimator is able to detect breaks near the boundary much more reliable than the likelihood-based approach. Furthermore, since we detect structural breaks by the Euclidean norm of the group of coefficient changes, the precision of the group LASSO estimator is mostly determined by the total magnitude of each break. This implies that we do not rely on a distinction between common and partial breaks which is important for the properties of hypothesis tests conducted in the likelihood-based approach to determine the number of breaks. In total, both approaches are conceptually very different so that one approach can serve as a valuable robustness check for the model specification chosen by the other approach.
Finally, we apply the two-step estimator to a term structure model of US interest rates to demonstrate its properties in a real world setting. Relying on the reduced computational burden of the proposed estimator, we are able to estimate the term structure model with daily data over a 30 year span and detect four important structural breaks. Our results reveal substantial differences in the parametrization of the resulting five term structure regimes. The estimates of the proportionality coefficient are smaller than one in every regime and the long-run implications of the expectations hypothesis are rejected in three out of five regimes.
The paper is organized as follows. (ref) outlines the proposed model selection procedure to estimate structural breaks in multivariate systems and presents our main technical results. (ref) is devoted to the Monte Carlo simulation study. (ref) reports the results of an empirical application of our methodology to the term structure of US interest rates, and (ref) concludes. Proofs of all theorems in the paper are provided in the Mathematical Appendix.
Using penalized regression techniques for structural break detection, we aim to divide a set of breakpoint candidates into two groups of active and inactive breakpoints. Our starting points are ChanYauZhang2014 and Schweikert2021, where a two-step procedure is proposed to detect and estimate multiple structural breaks in autoregressive processes and single equation cointegrating regressions, respectively. Here, the model of interest is a multiple equations system of linear regressions with integrated and stationary regressors, $q$ equations, and $T$ time periods.
We consider the following potentially cointegrated system in triangular form
where $Y_t$ is a $q \times 1$ vector of dependent variables, $X_t$ is a $r \times 1$ vector of integrated regressors, $w_t$ is a $s \times 1$ vector of stationary variables, $u_t$ and $\xi_t$ are I(0) error processes. The coefficient matrices $A$ and $B$ have dimension $q \times r$ and $q \times s$, respectively. Throughout, $\Vert \cdot \Vert$ represents either the Euclidean norm for vectors, i.e., $\Vert x \Vert = (\sum_{i=1}^n x_i^2)^{1/2}$ for $x \in \mathbb{R}^n$ or the Frobenius norm for matrices, i.e., $\Vert A \Vert = [\operatorname{tr} (A A')]^{1/2}$ for $A \in \mathbb{R}^{m \times n}$. We study the asymptotic properties of our estimator under the following assumptions about the involved processes:
The error term processes are assumed to be linear processes in Assumption (ref), satisfying the required conditions to ensure the validity of the functional central limit theorem for partial sum processes constructed from them (see, for example, Theorem 3.4 in PhillipsSolo1992 and its multivariate extension in Phillips1995). Further, $w_t$ is given as a stationary process with a sufficiently well-behaved distribution. Note that these assumptions could be replaced by other sufficient conditions that conform with a strong invariance principle or functional central limit theorem. We assume that $\Omega_{\xi}$ is positive definite which implies that $X_t$ is non-cointegrated. This ensures that the coefficients of $X_t$ have the standard $T$ rate of convergence. Assumption (ref) is chosen for technical reasons but it is less restrictive than it initially seems considering that $w_t$ might include the leads and lags of changes in $X_t$ (see, for example, Saikkonen1991, PhillipsLoretan1991, and StockWatson1993 for a treatment of second-order biased estimators in the context of cointegrating regressions).\footnote{We note that Assumption (ref) is stronger than the moment conditions given in Assumption (ref) (ii) and replaces them in the corresponding results (Theorem (ref) -- Theorem (ref)). However, Assumption (ref) is not necessary to proof Theorem (ref).}
Assumption (ref) allows for serially correlated error terms but this conflicts with endogenous regressors like lagged dependent variables. The results presented in the following remain valid for endogenous regressors as long as the errors are not permitted to be serially correlated BaiPerron1998.\footnote{The simulation results presented in Table S6 in Supplementary Material A show that the structural break detection remains consistent for endogenous regressors and white noise errors.}
We follow LiPerron2017 and OkaPerron2018 and apply scaling factors in Equation (ref) so that the order of all regressors is the same.\footnote{Higher order deterministic trends can be included in the model in the same way as long as the corresponding scaling factors are applied.} To simplify notation, we write the system of equations in its stacked form as
where $Z_t = (T^{-1/2} X_t', T^{-1}t, 1, w_t')'$ and $\theta = \operatorname{Vec}(A, \delta, \mu, B)$ is a $d = q(r + 2 + s)$ column vector, concatenating the coefficients for each regressor over all equations. The operator $\operatorname{Vec}(\cdot)$ stacks the rows of a matrix into a column vector, $\otimes$ denotes the Kronecker product and $I$ is a $q \times q$ identity matrix. While model (ref) allows for very flexible specifications and covers several special cases (e.g., SUR models for $A = 0$ and $\delta = 0$, VAR(s) models for $A = 0$ and $\delta = 0$ and $w_t = (Y_{t-1}, \dots, Y_{t-s})'$, or pure cointegrating regressions for $\delta = \mu = 0$ and $B = 0$), our main focus is on a full model specification with a mix of integrated and stationary regressors. For example, if $w_t$ contains the leads and lags of changes in $X_t$, we estimate structural breaks in a dynamically augmented cointegrating regression with multivariate responses. Assumption (ref) is violated in VAR models which invalidates our technical analysis for this specification. Hence, we refer to SafikhaniShojaie2020 for the appropriate assumptions to show that LASSO-type estimators can be used to estimate piece-wise stationary VAR models. If it is known to the researcher that some of the coefficients do not change, it is possible to introduce a selection matrix $S$ to consider only partial structural breaks in the sense that some coefficients are constant over the entire sampling period. This allows the researcher to estimate the respective coefficients with full efficiency. The selection matrix contains elements that are either 0 or 1 and hence, specifies which regressors appear in each equation.\footnote{Note that $S'S$ is idempotent with non zero elements only on the diagonal. The rank of $S$ is equal to the number of coefficients that are allowed to change.} For example, we can use the selection matrix to focus only on breaks in the matrix $A$, i.e., detecting breaks in the long-run coefficients of a cointegrating regression but leave the matrix $B$ constant over the sample period.
We assume that the system includes $m_0$ true structural breaks. Multiple (partial) structural breaks in the regression coefficients can be expressed using the following model
where $d(t_k) = 0$ for $t \leq t_k$ and $d(t_k) = 1$ for $t \geq t_k$. The total number of potential structural breaks or breakpoint candidates in this model is denoted by $m$, $\theta_{t_0}$ is the baseline coefficient vector, and $\theta_{t_k}$, $k = 1, \dots, m$ are regime-dependent changes in the regression coefficients. In situations where $m > m_0$, our structural break model in Equation (ref) considers more structural breaks than necessary which implies that the true coefficient vector does not change at some $t_k$, i.e., $\theta^0_{t_k} = 0$. For the breakpoints or break dates $t_k$, it holds by general convention that $1 = t_0 < t_1 < \dots < t_m < t_{m+1} = T + 1$. The relative timing of breakpoints is denoted by $\tau_k = t_k / T, k \in \lbrace 1, \dots m \rbrace$. We assume that there is a change in at least one coefficient matrix at each true structural break, so that $\Vert S \theta^0_{t_k} \Vert \neq 0$. To simplify the notation, we assume that all coefficients change at all breakpoints for the remainder of the paper.
In case of unknown number and timing of structural breaks, each point in time has to be considered as a potential breakpoint. Therefore, it is helpful to estimate the model in Equation (ref) with $m = T$ under the condition that the set $\theta(T) = \lbrace \theta_1, \theta_2, \dots, \theta_T \rbrace$ exhibits a certain sparse nature so that the total number of distinct vectors in the set equals the true number of breaks $m_0$. To use a convenient matrix notation, we define
$\mathcal{Y} = (Y_1', \dots, Y_T')'$, $\mathcal{U} = (u_1', \dots, u_T')'$ and $\boldsymbol{\theta}(T) = (\boldsymbol{\theta}_1', \dots, \boldsymbol{\theta}_T')'$. Furthermore, we define $\boldsymbol{Y} = \operatorname{Vec}(\mathcal{Y})$, $\boldsymbol{Z} = I \otimes \mathcal{Z}$, and $\boldsymbol{U} = \operatorname{Vec}(\mathcal{U})$. Now, the system for $m=T$ breakpoint candidates can be rewritten as
where $\boldsymbol{Y} \in \mathbb{R}^{Tq \times 1}$, $\boldsymbol{Z} \in \mathbb{R}^{Tq \times Td}$, $\boldsymbol{\theta}(T) \in \mathbb{R}^{Td \times 1}$, and $\boldsymbol{U} \in \mathbb{R}^{Tq \times 1}$. Note that this model specification reorders the groups so that $\boldsymbol{\theta}(T)$ contains the breakpoint candidates for each equation successively, i.e. $\boldsymbol{\theta}_i = \boldsymbol{K} \theta_i$ with commutation matrix $\boldsymbol{K}$. We assume that at least one of the baseline coefficients is nonzero to distinguish between active and inactive breakpoints without making any case-by-case considerations. This implies that the vector of true coefficients $\boldsymbol{\theta}^0(T)$ contains $m_0 + 1$ nonzero groups and $\boldsymbol{\theta}_1 \neq \boldsymbol{0}$. For the remainder of this paper, $\boldsymbol{\theta}_i = \boldsymbol{0}$ means that $\boldsymbol{\theta}_i$ has all entries equaling zero and $\boldsymbol{\theta}_i \neq \boldsymbol{0}$ means that $\boldsymbol{\theta}_i$ has at least one non-zero entry. We define the index sets $\bar{\mathcal{A}} = \lbrace 1 \leq i \leq T: \boldsymbol{\theta}^0_i \neq \boldsymbol{0} \rbrace$ denoting the indices of truly non-zero coefficients (including the baseline coefficient) and $\mathcal{A} = \lbrace i \geq 2: \boldsymbol{\theta}^0_i \neq \boldsymbol{0} \rbrace$ denoting the non-zero parameter changes. The set $\mathcal{A}_T = \lbrace \hat{t}_1, \hat{t}_2, \dots, \hat{t}_{\hat{m}} \rbrace$ denotes the $\hat{m}$ breakpoints estimated in the first step, i.e., indices of those coefficient changes which are estimated to be non-zero. $|\mathcal{A}|$ denotes the cardinality of the set $\mathcal{A}$ and $\mathcal{A}^c$ denotes the complementary set.
We propose to estimate the set of coefficient changes $\boldsymbol{\theta}(T)$ by minimizing the following penalized least squares objective function YuanLin2006:
where $\lambda_T$ is a tuning parameter and $\Vert \cdot \Vert$ denotes the $L_2$-norm. Minimizing the objective function in (ref) yields the group LASSO estimator which is denoted by $\hat{\boldsymbol{\theta}}(T)$. Using a group penalty for $\boldsymbol{\theta}_i$, in principle, we assume common breaks across all equations as defined by QuPerron2007 and OkaPerron2018. However, it should be noted that the issue of common or partial breaks is crucial for the specification of hypothesis tests required for the likelihood-based approach to detect the true number of breaks but it is not as important for the group LASSO estimator. Since we detect structural breaks by the Euclidean norm of the group of coefficient changes (vectorizing all coefficient matrices), only the total magnitude of each break enters our objective function. Consequently, prior knowledge about partial breaks would not improve our sensitivity detecting those breaks as much as it would in the likelihood-based approach.
Using these definitions, we frame the detection of structural breaks as a model selection problem and are able to use efficient algorithms from this strand of the literature HuangBrehenyMa2012, ChanYauZhang2014, YauHui2017 to eliminate irrelevant breakpoint candidates. Depending on the value of $\lambda_T$, a sparse solution is obtained so that the number of nonzero groups corresponds to the number of estimated breaks and the coefficient changes at each breakpoint are contained within the nonzero groups.\footnote{Practical guidance for the choice of $\lambda_T$ is given in (ref).} In the next subsection, we investigate the asymptotic properties of the first step estimator in this setting.
We show in the following that the group LASSO estimator for structural breaks in a system of linear regression equations is consistent in terms of prediction error but inherits the same problems, namely estimation inefficiency and model selection inconsistency, as shown for univariate AR models ChanYauZhang2014, single-equation cointegrating regressions Schweikert2021, and piecewise-stationary VAR models Gao2019. As discussed in ChanYauZhang2014, any two adjacent columns of the matrix $\mathcal{Z}$ only differ by one entry. Consequently, the restricted eigenvalue condition BickelRitovTsybakov2009 does not hold in our setting and we cannot establish our consistency proofs based on this assumption.
Further assumptions about the timing of true breakpoints ($\tau^0_k$, $k = 1, \dots, m_0$) and the magnitude of coefficient changes have to be stated to continue our analysis.
The first inequality of Assumption (ref)(i) is a necessary condition to ensure that a structural break occurs at $t_j^0$ and the second part excludes the possibility of infinitely large parameter changes.\footnote{Our definitions in Assumption (ref)(i) include the baseline coefficients which can, in principle, be relaxed but simplifies the technical analysis.} We do not consider breaks with local-to-zero behaviour in this setting (see BaiLumsdaineStock1998 for assumptions used in this context). This assumption is not believed to be restrictive for the intended empirical applications. In the case of a full specification with integrated regressors, applied researchers aim to estimate structural long-run relationships that are often needed for their follow-up analysis, e.g., in a cointegrated VAR as in Hansen2003. Essentially, they need optimal in-sample forecasts in terms of mean squared error of the cointegrating regressions under structural instability to consistently estimate deviations from the long-run equilibrium. BootPick2019 show that in-sample forecasts are largely unaffected by local-to-zero breaks. Assumption (ref)(ii) requires that the length of the regimes between breaks increases with the sample size albeit slower than $T$. If $\gamma_T$ is chosen with a slow enough rate, which depends on the tuning parameter sequence $\lambda_T$, it allows us to consistently detect and estimate the true break fractions.
The first result for the first step estimator shows that it is consistent in terms of prediction error if the tuning parameter $\lambda_T$ grows at the right rate.
The next result shows that the number of estimated breaks is at least as large as the true number of breaks. Furthermore, the location of the breakpoints can be estimated within an $T\gamma_T$-neighborhood of their true location. To state the theorem, we have to define the Hausdorff distance between the set of estimated breakpoints and the set of true breakpoints. We follow Boysen2009 and define $d_H(A, B) = \underset{b \in B}{\max} \, \underset{a \in A}{\min} |b - a|$ with $d_H(A, \emptyset) = d_H(\emptyset, B) = 1$, where $\emptyset$ is the empty set.
To obtain a consistent estimator for the number of breaks, their timing and coefficient changes, we need to design a second step refinement reducing the number of superfluous breaks. Immediate candidates are using a backward elimination algorithm (BEA) optimizing some information criterion ChanYauZhang2014, Gao2019, SafikhaniShojaie2020 or applying the adaptive group LASSO estimator as a second step using the group LASSO estimates as weights BehrendtSchweikert2020, Schweikert2021. In the following, we outline the former approach and discuss the latter approach in the Supplementary Material B because the BEA produces more accurate results in our simulation experiments.\footnote{Supplementary Material B can be found here: https://karstenschweikert.github.io/mequ_ci/mequ_ci_suppB.pdf}
According to Theorem (ref), the group LASSO estimator slightly overselects breaks under the right tuning. To distinguish between active and non-active breakpoints in the set $\mathcal{A}_T$, we employ an information criterion for the second step which consists of a goodness-of-fit measure, here the sum of squared residuals, and a penalty term as a function of the number of breaks. We define $\widehat{\widehat{\boldsymbol{\theta}}}_j$, $1 \leq j \leq m$ as the least squares estimator of $\boldsymbol{\theta}^0_j$, based on breakpoints estimated in the first step. Further, we define the sum of squared residuals over all $q$ equations as
where $\bar{Z}_t = (Z_t' \otimes I)$. For $m$ and the breakpoints $\boldsymbol{t} = (t_1, \dots, t_m)$, we can define the information criterion (IC)
where $\omega_T$ is the penalty term that is further characterized below in Theorem (ref). We estimate the number of breaks and the timing by solving
If the maximum number of breaks in the first step algorithm is chosen to be small, the evaluation of the information criterion for each combination of breakpoints can be achieved easily. The following result shows that minimizing the IC gives us a consistent estimator for $m_0$ and $\mathcal{A}$.
If $|\mathcal{A}_T|$ is relatively large, evaluating every combination of breakpoints again becomes computationally intensive. Hence, we follow ChanYauZhang2014 and use a backward elimination algorithm to successively remove the most redundant breakpoint that corresponds to the largest reduction of the IC until no further improvement is possible. This means that including this breakpoint improves the fit sufficiently to outweigh the costs of estimating the coefficients for an additional regime. The details of the algorithm are outlined in Supplementary Material A.\footnote{Supplementary Material A can be found here: https://karstenschweikert.github.io/mequ_ci/mequ_ci_suppA.pdf} We denote the set of estimated breakpoints obtained from the BEA by $\mathcal{A}^*_T = (\widehat{t}^*_1, \dots, \widehat{t}^*_{|\mathcal{A}^*_T|})$. The next theorem shows that the estimator based on the BEA has identical asymptotic properties.
Using the BEA, it is also possible to optimize another information criterion, say the BIC, for each regime to eliminate some variables from all equations. This allows us to investigate whether some variables lose importance during parts of the sample period. Depending on the chosen model structure this could even be interpreted as some variables dropping out of the long-run equilibrium relationship for a certain period.
Applying scaling factors to the integrated regressors and the linear trend in our model ensures that all regressors have the same order. In turn, this means that the OLS estimator in the second step, $\widehat{\widehat{\boldsymbol{\theta}}}_j$, has the same convergence rates for all coefficients in the model. In principle, it is also possible to conduct post-LASSO OLS estimation (after the second step) without scaling factors to benefit from higher convergence rates of the estimator for coefficients of trending variables.
We conduct simulation experiments to assess the adequacy of our technical results presented in (ref). Specifically, we investigate the finite sample performance of our estimator with respect to the speed and accuracy in finding the exact number of breaks and their location. In principle, we can optimize (ref) directly using a coordinate descent algorithm like the one proposed in Breheny2015 for a grid of $\lambda_T$ values. However, such an algorithm is not computationally efficient for change-point problems. Thus, for the simulations and the empirical application, we rely on a modified group LARS algorithm that approximates the solution for the first step and the BEA for the second step. The choice of the tuning parameter $\lambda_T$ is translated to pre-specifying the maximum number of breakpoint candidates $M$, i.e.\ the maximum number of non-zero groups in $\hat{\boldsymbol{\theta}}(T)$ that is returned when the LARS implementation is used. The modified group LARS algorithm evaluates each point in time as a breakpoint candidate but returns only the $M$ most relevant breaks. According to Theorem 2, we know that the group LASSO estimator overselects breaks in the first step. Hence, $M$ should be set large enough to encompass all true breakpoints and some additional falsely selected non-zero groups. The BEA then asymptotically guarantees that the set of change-points is attained in the second step. Further, the minimum distance between breaks needs to be specified based on the number of coefficients in the model to guarantee that these coefficients can be estimated accurately in each regime. Details about the modified group LARS algorithm, the BEA, and additional simulation results are included in Supplementary Material A.
We consider model specifications with one, two and four breakpoints, respectively. The following DGP is employed to model a multiple equations cointegrating regression with multiple structural breaks,
where $X_t = (X_{1t}, X_{2t}, \dots, X_{Nt})'$, $\Sigma_{u,ii} = \sigma^2_u$, $\Sigma_{\xi,ii} = \sigma^2_{\xi}$ and $\Sigma_{e,ii} = \sigma^2_e$, for $i = 1, \dots, q$, i.e.\ the innovations of our generated processes have multivariate normal distributions with identical variances. We choose the variances of these processes so that each column of $Z_t = (T^{-1/2} X_t', T^{-1}t, w_t')'$ has variances less than or equal to the error term variances. To achieve this, we set the variances to $\sigma^2_{\xi} = \sigma^2_e = \sigma^2_u = 1$. $\mu$ is a non-zero intercept vector, $A_t$ and $B_t$ are time-varying coefficient matrices with at least one non-zero entry in the baseline specification and a finite number of breaks. $\delta_t$ is a time-varying $q$-dimensional vector and $\Phi$ is a coefficient matrix for the VAR(1) process that fulfills the required stationarity conditions. For simplicity, we choose a diagonal matrix with each diagonal entry equal to 0.5. For the main results, we set $q = 2$ and use the following coefficient matrices:
with $c = 1$ and the subscript $t_0$ denoting the initial coefficient matrix before the first breakpoint. Moreover, we set $\mu = (2,2)'$, $\delta_{t_0} = (2,2)'$ and $\delta_{t_j} = \delta_{t_j-1} + c (2,2)'$ for $j = 1, \dots, m_0$. In our baseline specification, the coefficient changes amount to two standard deviations of the error terms. Similar to BaiLumsdaineStock1998, we also specify $c \in \lbrace 0.25, 0.5, 1.5 \rbrace$ to investigate the performance for smaller and larger break magnitudes. Moreover, we investigate the effect of cross-correlated errors, serially correlated errors, and endogenous regressors. The results of those robustness checks are included in Supplementary Material A to conserve space.
Naturally, the ability of all structural break estimators to detect breaks depends on the overall signal strength. NiuHaoZhang2015 define signal strength in change-point models by $S_{NHZ} = m_{\theta}^2 I_{\min}$, where $I_{\min} = \min_{1 \leq j \leq m_0+1} \vert t_j - t_{j-1} \vert$ is the minimum distance between breaks and $m_{\theta}$ is the minimum jump size as defined in Assumption (ref). For our main simulations concerned with showing the consistency of the two-step estimator, we follow Schweikert2021 and use equal jump sizes for multiple breaks as well as locating the breaks with equidistant spacing between them. Hence, overall signal strength is a linear function of the sample size in our simulations. We choose a minimum of 50 observation per regime and double the sample sizes in line with the conventional asymptotics specified in (ref). Consequently, the sample sizes chosen differ for an increasing number of breaks. Note that 12 coefficients are present in the full model specification which requires a substantial number of observations in each regime to estimate them precisely.
In (ref), we report our results for $r = 2$ integrated regressors, $s = 2$ stationary regressors, a time trend, and $q = 2$ equations. We specify our model for one break located at $\tau = 0.5$, two breaks at $\tau = (0.33, 0.67)$ and four breaks at $\tau = (0.2, 0.4, 0.6, 0.8)$ to have an equidistant spacing on the unit interval. Since the choice of a change-point detection algorithm is often reflective of a trade-off between speed and accuracy, we first compute the average computing time in seconds (Time [s]) for each sample size.\footnote{All simulations experiments are conducted on a computer with an Intel i5-6500 CPU at 3.20GHz and 16GB RAM.} Next, we compute the percentages of correct estimation (pce) of the true number of breaks $m_0$ and measure the accuracy of the break date estimation conditional on the correct estimation of $m_0$. For this matter, we compute the standard deviations of the estimated relative timing. The estimated coefficients are not reported to conserve space but can be obtained from the author upon request.
Our simulation results reveal that the two-step estimator (panel A) is less precise compared with the likelihood-based approach (panel B) when it comes to the estimation of the individual break locations. This is not surprising considering that the likelihood-based approach according to QuPerron2007 rests on a dynamic programming algorithm and uses repeated OLS regressions to determine the location of the breaks in an almost exact fashion. The trade-off here is then computing time versus precision. As outlined above, the computational efficiency of the two-step estimator is much higher ($O(M^3 +MT)$ versus $O(T^2)$) which manifests in reduced computing times at moderate to large samples. For example, while it takes the likelihood-approach several minutes (on average 2,013 seconds) on a modern computer to solve the change-point problem for $T=2,000$, the two-step estimator solves the same problem in seconds (on average 36 seconds). It is also important to note that, in contrast to the two-step estimator, the likelihood-based approach is conceptually not able to consistently estimate the exact number of breaks. As expected, the sequential test procedure reaches its specified nominal confidence level at large samples sizes and in our setting ($\alpha = 0.05$) detects the number of breaks in roughly 95% of all cases.\footnote{It seems to be undersized for small sample sizes and a larger number of breaks, reaching a 100% detection rate for $m = 4$ and $T = 250$ but then declining to 95% for larger samples.} Instead, the two-step estimator already attains close to a 100% detection rate at moderate sample sizes. This also means that the most difficult cases (leading to a non-rejection of the sequential test's hypothesis) are not considered for the evaluation of the precision of the likelihood-based approach because in columns two to six of (ref) only those cases with the correctly estimated number of breaks can be properly evaluated. BaiPerron1998 suggest to let the size of those sequential tests go to zero asymptotically to avoid misspecification. However, to implement this in practice is difficult without any guidance on the exact rate for an adjusted $\alpha$. Alternatively, an information criterion might be used to compare specifications with a different number of breaks.
We further study the performance of the two-step estimator in settings with structural breaks of smaller or larger magnitude.\footnote{The results can be found in Tables S1 -- S3 in Supplementary Material A.} For instance, we apply the factor $c=1.5$ to increase the break magnitude in Equation (ref) to three standard deviations of the error term. Although we find a slightly better detection of the number of breaks, we observe no substantial effect on the precision in terms of finding the true location of the breaks. It appears that increasing the magnitude does not further improve the performance in small samples. A reason for this upper bound in precision is the pre-selection of breakpoint candidates in the first step which is only accurate up to a $T\gamma_T$ neighborhood of the true breakpoints. In contrast, reducing the magnitude when applying the factor $c=0.5$, i.e., to one standard deviation, leads to worse results for small sample sizes. For example, in the case of four active break and $T=250$ observations, we find that the rate of correctly detecting the number of breaks drops from 89.0% in the baseline specification to 68.3%. However, when the sample size increases to $T=500$, we almost reach the same rate of detection and a similar precision for the break's location. The likelihood-based approach is still very precise when the break magnitudes are reduced by the factor $c=0.5$. However, when we reduce them further by setting $c=0.25$, the detection rate breaks down to, e.g., 6.7% for four breaks and $T=250$ observations while the two-step approach still reaches a 40.1% detection rate. We conclude that the performance of the two-step approach naturally depends on the size of the break magnitude but it is fairly robust in medium to large sample sizes which should be its primary field of application considering its improvements in terms of computational costs.
To investigate whether the results are driven by the integrated regressors or the linear trend in the model, we run the simulation experiments for a reduced SUR model specification. Here, we generate data under the restriction that $A_i = diag(0)$ and $\delta_i = 0$ for $i = 0, \dots m_0$. The results are reported in (ref). We consider two variants of the simulation experiment: first with no cross-correlation between the error terms for both equations as in Equation (ref) and second with cross-correlation coefficient $\rho = 0.5$. We find almost identical precision for the SUR model and no substantial difference in the presence of moderate cross-correlation. This shows that the precision of the two-step approach is predominately driven by the magnitude of breaks in terms of their Euclidean distance for the vector of coefficients and the magnitude is the same for the SUR specification. The reduced number of coefficients in the SUR model does not seem to improve the performance substantially, because the sample size in each regime is already large enough to estimate the full model. Additional simulation experiments (not reported) show that this aspect certainly gains importance for smaller regimes with less observations.
We extend the model with a third equation to investigate if the results either improve because a common break is indicated in another regression equation BaiLumsdaineStock1998, QuPerron2007 or deteriorate because the detection relies on additional coefficient changes that need to be estimated. We consider a special case in which the break magnitudes stay constant after adding the third equation. In practice, we hope that estimating the structural breaks jointly in all equations includes some larger coefficient changes that help to find the common break dates. The results for the $q=3$ case, reported in (ref), show that we reach similar detection rates and approximately the same precision for the timing of the breaks. As the break magnitudes are identical for each coefficient, the $q=3$ setting is a straightforward extension of the $q=2$ setting. To investigate further whether breaks in a subset of the coefficients can be reliably detected in a larger system, we consider a partial break specification in the $q=2$ case. Here, only the coefficients of the first equation change. The results are reported in Table S5 in Supplementary Material A and again show that the two-step estimator's performance is mostly determined by the total break magnitude which is slightly reduced in this case.
Finally, we turn to the one-break case of our main specification again and move the break to the right boundary of the unit interval. We fix the number of observations in the second regime to 25 and increase the overall sample size. Consequently the breakpoint moves further to the boundary each time we double the sample size. The true break fractions for the sample sizes $T \in \lbrace 100, 200, 400, 800 \rbrace$ are 0.75, 0.875, 0.9375, and 0.96875, respectively. We also consider a two-break setting where the middle regime shrinks relative to the total sample size. We set the first breakpoint at $\tau_1 = 0.5$ and the second one at $\tau_2 = 0.5 + 25/T$. The results for both experiments are reported in (ref). These exercises show that the two-step estimator is robust to setting a small minimum break distance while the sequential tests to detect breaks in the likelihood-based approach require a rather large trimming parameter to ensure the right size and sufficient power. Of course, if the trimming parameter is set to 0.15 as it is suggested for empirical data QuPerron2007, it is infeasible to detect breaks that are too close to each other or near the boundary of the unit interval. In case of one break at the boundary, the likelihood-based estimator in our simulations (almost) always falsely indicates a break at 0.85.
In our empirical application, we apply the two-step estimator to US term structure data. Thereby, we revisit the study by Hansen2003 who proposes a framework to test for structural change in cointegrated vector autoregressions. The author tests two potential structural breaks in September 1979 and October 1982 that coincide with large changes in the Fed's policy. Only after accounting for these structural changes, the long-run implication of the expectations hypothesis (EHT) cannot be rejected. In the subsequent analysis, we study more recent data on the US term structure and detect structural breaks without assuming any prior knowledge about their location.
Following Campbell1987, we expect that the term structure of interest rates is determined by the expected future spot rates being equal to the future rate plus a time-invariant term premium in the long-run. This implies that, independent of the maturity, the yields should be cointegrated with pairwise cointegrating vector $(1,-1)$. However, several early studies report that the EHT fails in empirical practice Froot1989, Campbell1991. Besides other empirical difficulties, structural breaks are named as one of the important reasons for this failure Lanne1999, Sarno2007, Bulkley2011. Several studies investigate whether regime shifts in the term structure of interest rates are related to changes in monetary policy Tillmann2007, Thornton2018. The important question for applied researchers and policy-makers is whether these equilibrium relationships are robust over different regimes.
Using daily data from January 1990 to July 2021 on the term structure of US interest rates, we end up with more than 8,000 observations to estimate the term structure model. We again emphasize that almost exact segmentation algorithms like the one used for the likelihood-based approach are substantially slower than the two-step procedure based on the group LASSO estimator (26,380 seconds versus 25 seconds for solving the change point problem). Taking into account that several re-estimations of the model must be conducted to find the right specification and perform robustness checks, a reduced computational burden is important to encourage routine checks for structural breaks in multivariate systems. We use fitted yields on zero coupon US bonds with 10-year ($r_{10y,t}$), 5-year ($r_{5y,t})$, and 1-year ($r_{1y,t}$) maturity in the term structure model,
We choose a model specification matching the long-run component in the cointegrated VAR used in Hansen2003. Each equation models the pairwise relationship between the longer term maturity and the short term maturity (a one-to-one relationship under the EHT), while the constant accounts for a term premium. Additional maturities could be analyzed, leading to additional equations in the model, but we try to maintain a simple model structure. The data are produced according to the approach of Kim2005 fitting a simple three-factor arbitrage-free term structure model to U.S. Treasury yields since 1990, in order to evaluate the behavior of long-term yields, distant-horizon forward rates, and term premiums.\footnote{The data can be downloaded from the St Louis Fed's database FRED.} (ref) provides a time series plot of the data. We assume that the individual variables follow unit root processes.\footnote{See the discussion in Hansen2003 and the references given therein why this assumption is useful for the empirical modelling of the term structure although some features like the non-negativity of interest rates are arguments against it. We also conduct unit root tests which do not provide evidence against this hypothesis in this specific sample period.} The results of a Johansen trace test for the full sample suggest that the trivariate system is cointegrated with cointegration rank one or two depending on the specification of the deterministic terms. According to the long-run implications of the EHT, we would expect a cointegration rank of two but it is well-known that the Johansen trace tests are not robust to structural breaks in the cointegrating vectors and the deterministic terms Lutkepohl2004, Saikkonen2006.
We estimate the cointegrated term structure regression in Equation (ref) with a dynamic OLS specification adding two leads and lags of $\Delta r_{1y,t}$. First, we estimate the model for the full sample without accounting for any structural breaks. The coefficient estimates are $\hat{\beta}_1 = 0.764 (0.116)$ and $\hat{\beta}_2 = 0.894 (0.078)$ which lead to a rejection of the EHT at the 5% significance level for the 10-year maturity but not for the 5-year maturity.\footnote{Bootstrap standard errors are computed based on 600 replications of the sieve bootstrap method for cointegrating regressions proposed in Chang2006.} In a second step, we try to capture all relevant structural breaks. Due to the large number of observations, we pre-specify a large maximum number of breaks, $M=40$, and maintain a minimum break distance of two month (50 daily observations) to obtain accurate coefficient estimates in each regime.\footnote{The results are robust for different choices of $M$ as long as $M > 4$. Since a larger value of $M$ allows the modified group LARS algorithm to take more steps, a larger value of $M$ can, in principle, positively affect the precision of the estimated break locations. However, the estimates do not change for $M > 40$.} We consider a specification with a constant and apply the required scaling factors so that each regressor has the same order.\footnote{The results with and without linear trend do not differ substantially.} Using the two-step group LASSO estimator, we obtain four structural breaks. All breakpoint estimates can be related to important monetary policy events. To show that, we depict the trajectory of the effective federal funds rate (EFFR) in (ref) and indicate the regimes. The first break is located in November 1994 after the EFFR sharply increases and the yield spread narrows, the second break is located in April 2003 during the recession. Here, the EFFR falls dramatically and we observe wider spreads. The third break is located in August 2010 after the Global Financial Crisis, and the fourth break is located in March 2015 at the end of the zero target rate regime. After the structural breaks are obtained, we re-estimate the model in each regime without scaling factors for higher precision and report the resulting coefficients in (ref). For comparison, the estimated breakpoints from the likelihood-based approach are located in January 1995, December 2003, July 2008, and July 2013 corresponding to roughly the same monetary policy events. The algorithm also detects four breakpoints with the sequential test clearly rejecting the hypothesis of three breaks in favor of a four break model. More than four breaks cannot be allocated if the default minimum regime length should be maintained (setting the trimming parameter to 0.15).
It can be observed that the pairwise cointegrating vectors for most regimes substantially differ from $(1,-1)$. A simple $t$-test of the hypothesis does not lead to a rejection of the EHT for the last two regimes.\footnote{Note that the reported standard errors do not take the break estimation uncertainty into account.} While the estimates of the proportionality coefficient for the first three regimes in case of the 10-year maturity are reasonable (albeit still being able to reject the EHT), the estimated coefficient for the fourth regime from August 2010 to March 2015 is negative with a much larger standard error. The non-rejection of the EHT in this regime can therefore be attributed to the higher error term variance in this period. This regime can also be associated with a zero target rate and an unusually steep yield curve which might explain the unusually small proportionality coefficient for the 10-year maturity. The final regime is marked by relatively narrow yield spreads that begin to widen during the COVID-19 pandemic. Still, our results suggest that this most recent regime comes closest to satisfying the EHT. In summary, accounting for multiple structural breaks in the term structure model reveals some important differences for subsamples of the data but does not solve the term structure puzzle. Instead, it seems that the EHT does not hold for most of the sampling period.
We have proposed a computationally efficient alternative to the existing likelihood-based approach solving the change-point problem in multivariate systems with a mix of integrated and stationary regressors. Our two-step estimator is able to consistently determine the number of breaks, their timing and to jointly estimate the coefficients for each regime. This can be achieved without the need to conduct sequential tests and does not require generating critical values beyond those supplied in the original paper, for example, in situations with larger systems composed of many regressors. The algorithm is much faster than the dynamic programming algorithm used in the likelihood-approach and can solve change-point problems for large multivariate systems and several thousand observations in seconds. In turn, the likelihood-based approach allows for a straightforward construction of valid confidence bands by inverting the likelihood ratio test for a given break date EoMorley2015. It remains to be investigated if such confidence bands can be constructed within the model selection approach taken in this paper.
The crucial first step estimation is based on the group LASSO estimator and we utilize a group LARS algorithm to solve the change-point problem. Alternative choices of other penalties in the objective function could be used to potentially improve the estimator in some directions. For example, SafikhaniShojaie2020 use a fused LASSO penalty to allow the number of equations to grow with the sample size. Such high-dimensional extensions might be useful in panel settings where both the time and the cross-sectional dimension increase asymptotically. In principle, the proposed two-step estimator can also be applied to change-point problems in other important model classes. If we relaxed the strict exogeneity assumptions and instead used more restrictive assumptions about the error terms, for example, assuming Gaussian white noise errors, we could also deal with structural breaks in VECMs which involve a mix of integrated and stationary variables. Since this application poses additional technical difficulties, we leave this topic for future research.
The author thanks Konstantin Kuck, Thomas Dimpfl, Markus M\"o\ssler, Robert Jung, and participants of the Research Seminar in Economics at the University of Hohenheim, the Asian Meeting of the Econometric Society in China 2022, the Econometric Society 2022 Australasia Meeting, the AMES 2022 Tokyo, and the EEA-ESEM in Milan for valuable comments. Further, he thanks Chun Yip Yau for sharing the programs for the group LARS algorithm, and Zhongjun Qu, Pierre Perron and Tatsushi Oka for sharing the programs for the likelihood-based approach. Funding by the German Research Foundation (Grant SCHW 2062/1-1) is gratefully acknowledged.