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.
50,415 characters · 11 sections · 69 citation commands
Change-Point Detection in Time Series Using Mixed Integer Programming
\def\spacingset#1{ {#1}} \spacingset{1} \if00 {
} \fi \if10 \fi
{\it Keywords:} Structural breaks, $l_0$-penalization, $l_1$-penalization, mixed integer quadratic programming \spacingset{1.8}
Identifying structural breaks in time series, also known as change-points, regime shifts and concept drifts, is a major area of interest within theoretical and applied statistics, going back at least to the 1960s Shiryaev, Roberts. In modern econometrics the focus has been on statistical approaches that estimate breakpoints by minimizing the regression sum of squares bai1998estimating, bai2003computation or $l_1$-penalized sum of squares qian_su_2016, kaddoura/westerlund:23. As one of the most highly-cited examples, bai1998estimating proposed a specific-to-general testing strategy for estimating the number of breaks in linear regression models with potential heterogeneity in the errors. The method requires testing the null hypotheses of $m$ breaks against the alternative of $m+1$ break starting with $m=0$. The estimated number of breaks then is that for which the null hypothesis is not rejected.
Testing-based approaches have been criticized for not always offering a consistent estimator of break dates and for the tendency to overestimate the true number of breaks with a positive probability, equal to the tests' significance level asymptotically. bai2003computation suggested using information criteria to choose the number of breaks, providing a consistent estimator of the break number. However, the approaches inevitably put restrictive assumptions on the minimal length of a regime to be set by the researcher, while the critical values crucially depend on this length both in large and in small samples.
Penalized methods have been proposed to circumvent the restriction on the minimal length of regimes. The LASSO (Least Absolute Shrinkage and Selection Operator) of tibshirani1996regression has been extremely effective at selecting the number of regression parameters with a simultaneous estimation of the non-zero parameters in linear models. harchaoui2010multiple and bleakley2011group consider the estimation of break locations in one-dimensional piece-wise constant signals, under the assumption of independence. chan2014group extend their approach to dependent data allowing the number of breaks to grow with the sample size. They also provide a justification for using a second step in the selection procedure in order to prevent an overestimation of the number of breaks.
behrendt2021note proposed using adaptive group LASSO to select the number of breaks consistently as an alternative to the two-step procedure of chan2014group. The two-step procedure of chan2014group is easier to apply, but it is less efficient than adaptive group LASSO of behrendt2021note. qian_su_2016 also considered a linear regression model and estimated the number of regimes and model parameters by using adaptive fused LASSO. Their approach is also two-step due to an overestimation of the true break date in the first step.
In the context of these developments, it has been a common belief that mixed integer optimization (MIO) is not suitable for such problems due to what is known as “combinatorial explosion”, that is, the explosive growth in the number of combinations to consider and associated insurmountable computational task. However, recent remarkable advances in computational and algorithmic methods of optimization over integer-valued arguments have shown attractive properties of integer and mixed integer programming as a means of obtaining efficient and provably optimal solutions in a wide range of statistical problems bertsimas2016, mazumder2023subset, hazimeh2023grouped, gomez2021mixed, rebennack2020piecewise.
For example, the problem of subset selection has become feasible and even standard in applications with datasets much larger than statisticians previously thought possible bertsimas/etal:20. Moreover, it is often suggested that subset selection using integer programming outperforms LASSO in many situations. The MIO challenge was famously picked by statisticians in a recent issue of Statistical Science hastie/etal:20. Such developments raise the prospect of applying MIO in other settings of interest to econometricians. For example, lee/etal:17 propose and implement MIO for a change-point regression model with one change and apply it to model the US unemployment rate.\footnote{We thank an anonymous referee for pointing out this paper to us. The paper considers a similar model to ours, but with only two regimes and the difference in the slope parameters determined by possibly unobserved factors. It also designs advanced computational methods for when the dimension of $\beta$ grows. In contrast, our approach focuses on the estimation of constant slope coefficients with an unrestricted number and location of breaks.}
In this paper we develop an MIO-based framework for simultaneously estimating the number and location of structural breaks as well as the parameters of a time series regression model. Formulated as a mixed integer quadratic programming problem, the new method shows an attractive performance compared to LASSO-type procedures, especially when the number of breaks is high. Existing solvers quickly obtain solutions of the problem to optimality, and no multi-step procedures, typical for LASSO-based competitors, are necessary. Importantly, our framework permits estimation of the unknown number of structural breaks while accommodating the option of specifying a required or minimal number of breaks if necessary. As we show, the new estimator enjoys relative simplicity and attractive theoretical properties under assumptions that are somewhat weaker than in the alternatives available in the literature.
To demonstrate the effectiveness of our approach, we conduct comprehensive Monte Carlo simulations, using a well-known and easily accessible solver. We show how to choose the tuning parameters and compare our results with established methods such as those proposed by bai1998estimating and qian_su_2016.
The paper is organized as follows. Section (ref) formulates the model and assumptions, and proposes our MIO approach. In Section (ref), we state the additional assumptions and main asymptotic results on consistency and asymptotic normality of the proposed estimators. Numerical experiments are discussed in Section (ref). Empirical applications are provided in Section (ref). Section (ref) concludes. All proofs and technical details are collected in the Appendix.
We assume that the data is generated by the following process:
where $x_t$ is a $p\times 1$ vector of regressors, $u_t$ is the error term, and the $(p\times 1)$ vector $\boldsymbol \beta^*_t$ takes distinct vector values $\boldsymbol \alpha^*_j, j=1, \ldots, m^*+1$, in the time interval $T^*_{j-1} \le t < T^*_{j}$, where $m^*$ is the number of breakpoints and where we use the convention that $T^*_{0}= 1$ and $T^*_{m^*+1}= T+1$. In this model, the indices $(T^*_1,\dots,T^*_{m^*})$, or breakpoints, are assumed to be unknown.
The goal is to find the unknown number $m^*$ of unknown break dates $(T^*_1,\dots,T^*_{m^*})$ as well as the regression coefficients $\boldsymbol \alpha=({\boldsymbol \alpha^*_1}^{\top},\dots, {\boldsymbol \alpha^*_{m^*+1}}^{\top})^{\top}$. It is clear that with no penalty, an in-sample prediction error minimization for ((ref)) gives $\hat{m}=T-1$ breaks and a perfect fit. This solution is unlikely to generalize well out-of-sample. To avoid overfitting, it is natural to impose a penalty that counteracts the reduction in prediction error for adjacent values of $\boldsymbol \beta^*_t$ that are not too far from one another. A common way of doing this is to utilize various forms of $l_1$-norm of the difference $\boldsymbol \beta^*_t-\boldsymbol \beta^*_{t-1}$; see, e.g., Group Fused LASSO (GFL) of qian_su_2016, Grouped LASSO of kaddoura/westerlund:23.
The class of estimators we consider can be stated as the following $l_0$-penalized optimization problem:
where $\mathbf{1}\{\cdot\}$ denotes the indicator function and $\lambda=\lambda_T$ is a positive tuning parameter. We note that our penalty function, $\sum_{t=2}^{T}{\mathbf{1}\{\boldsymbol \beta_t\ne \boldsymbol \beta_{t-1}\}}$, simply counts the number of jumps in the regression coefficient vectors $\boldsymbol \beta_t$.
It is interesting to compare and contrast our approach with the GFL method of qian_su_2016, which is based on the following optimization problem:
The difference between the two estimators is in the penalty function. While our approach penalizes the number of jumps in the $\boldsymbol \beta$-coefficients directly, the GFL method does so indirectly, by using a group LASSO-type penalty $\sum_{t=2}^{T}\|\boldsymbol \beta_t- \boldsymbol \beta_{t-1}\|$. This penalty serves as a proxy for our $\ell_0$ penalty $\sum_{t=2}^{T}{\mathbf{1}\{\boldsymbol \beta_t\ne \boldsymbol \beta_{t-1}\}}$, which counts the exact number of coefficient jumps. We note that one of the consequences of using a LASSO-type ($\ell_1$) penalty rather than an $\ell_0$ penalty is the resulting shrinkage of the estimated coefficient differences $\boldsymbol \beta_t- \boldsymbol \beta_{t-1}$.
To leverage impressive recent advances in the field of mixed integer optimization (MIO), we propose solving ((ref)) by formulating it as a MIO problem as follows:
where $\boldsymbol{e}$ is a $p$-vector of ones. Appendix A contains the details of why the two formulations are equivalent. Formulation ((ref))-((ref)) involves optimization over continuous variables $\boldsymbol \beta$ and binary variables $z_t$, where $z_t$ equals 1 if and only if there is a break at time $t$. In ((ref)), $M$ is some very large pre-specified constant, which gives name to the “Big-M” term for the general formulation of the above problem. We note that Big-M formulations are extensively used in mixed integer programming. In particular, they have recently been used in a number of $\ell_0$-regularized regression approaches bertsimas2016,mazumder/radchenko:17,hazimeh2023grouped. Optimization problem ((ref))-((ref)) is tractable and can be efficiently solved in practice using state-of-the-art MIO solvers.
In this section we study the asymptotic properties of our proposed estimator ((ref)). We start by introducing some notation. Let $I^*_j=T^*_j-T^*_{j-1}$ for $j=1,...,m^*+1$ and define $$ I_{\text{min}}=\min_{1\le j\le m^*+1}I^*_j, \quad J_{\text{min}}=\min_{1\le j\le m^*}\|\boldsymbol \alpha^*_{j+1}-\boldsymbol \alpha^*_{j}\|, \quad \text{and} \quad J_{\text{max}}=\max_{1\le j\le m^*}\|\boldsymbol \alpha^*_{j+1}-\boldsymbol \alpha^*_{j}\|. $$ We note that $I_{\text{min}}$ is the smallest interval length among the $m^*+1$ regimes of the true data-generating process, while $J_{\text{min}}$ and $J_{\text{max}}$ measure the smallest and largest jump sizes, respectively, in the true vector of coefficients.
The main result in this section establishes consistency of our approach in estimating the true number of breaks, breakpoints, and regression coefficients and also derives the corresponding rates of convergence. This result corresponds to the combination of the following two theorems in qian_su_2016: Theorem 3.4 (on correctly estimating the true number of breaks using the information criterion) and Theorem 3.1 (on the rate of convergence for the breakpoints and coefficients when the correct number of breaks is used). We impose the same assumptions (A1 and A2) on the $\{(x_t,u_t)\}$ process as qian_su_2016 do in their theoretical analysis. In particular, we assume that $\{(x_t,u_t)\}$ is a strong mixing process with a geometric decay rate and impose moment conditions on the corresponding random variables. We also bound the eigenvalues of $\frac1{r-s}\sum_{t=s}^{r-1} E(x_tx_t^{\top})$ for $r-s\ge T\delta_T$ and lower-bound $\delta_T$ depending on whether the finite-moment or the exponential-moment condition is satisfied. More specifically, if $(x_t,u_t)$ have uniformly bounded $4q$-th moments for some $q>1$, then $\delta_T$ cannot go to zero faster than $T^{(1-q)/q}$; alternatively, if $Ee^{c(\|x_t\|+|u_t|)^{2\gamma}}$ are uniformly bounded for some positive $c,\gamma$, then the fastest rate of decrease for $\delta_T$ is $(\log T)^{(2+\gamma)/\gamma}/T$. Assumptions A1 and A3 are formally stated in Appendix B. We also impose the following additional requirements.
\noindentAssumption A3.
This assumption is a weaker version of Assumption A3 in qian_su_2016, which is required for their Theorem 3.4. More specifically, we do not impose their conditions $Tm^*\big[(\log I_{\textrm{min}})^{c_{\delta}/2}T^{-1/2}I_{\textrm{min}}^{-1/2}+I_{\textrm{min}}^{-1} \big]\big( I_{\textrm{min}} J^2_{\textrm{min}} \big)^{-1}\rightarrow0$ and $m^*=O\big(\log T\big)$.
We define $\widehat{m}$, $\widehat{\boldsymbol \alpha}_j$, $\widehat{T_j}$ as the characteristics of our estimator ((ref)) that are the analogs of the corresponding population quantities ${m}^*$, ${\boldsymbol \alpha}^*_j$, and ${T}^*_j$.
\noindentProof: See Appendix C.
In the two remarks above, we compare our estimator to the following two-stage procedure of qian_su_2016. First, a base GFL estimator is obtained for a range of values of the tuning parameter $\rho_T$; second, the final estimator is determined by selecting the tuning parameter using an information criterion that penalizes the number of breaks. In contrast, our approach does everything in one go, and avoids the estimation bias that comes from LASSO penalty.
We can see that, in comparison to the base GFL estimator of qian_su_2016 our estimator enjoys better asymptotic properties. In particular, while qian_su_2016 show that the GFL estimator has at least as many breaks as the true model (see their Theorem 3.3), they do not establish a complimentary upper bound result. In contrast, we show that our estimator recovers the correct number of breaks with probability tending to one.
Next, we establish the asymptotic normality of our estimated regression coefficients. To state a clean result, we assume that $m^*$ is fixed and nonzero. However, we note that this result can be extended to the general case as in qian_su_2016, by imposing additional assumptions on $m^*$ and stating the central limit theorem for pre-specified fixed-dimensional sub-vectors of coefficients.
We impose the following additional conditions, which are also required by qian_su_2016 in the analogous result for their estimator.
\noindentAssumption A4.
Because our estimator recovers the correct number of breakpoints with probability tending to one, we follow the approach of qian_su_2016 and establish asymptotic normality for the estimator that solves optimization problem (2) with the restriction that the total number of breakpoints is exactly $m^*$, i.e., $\sum_{t=1}^T z_t=m^*$. We write $\widehat{\boldsymbol \alpha}_{m^*}$ for the corresponding vector of estimated regression coefficients and observe that $\widehat{\boldsymbol \alpha}_{m^*}=\big(\widehat{\mathbb{X}}^{\top}\widehat{\mathbb{X}}\big)^{-1}\widehat{\mathbb{X}}^{\top}Y$, where $\widehat{\mathbb{X}}=\text{diag}\big( (x_1,...,x_{\widehat{T}_1-1})^{\top},...,(x_{\widehat{T}_{m^*}},...,x_{T})^{\top} \big)$ and $Y=(y_1,...,y_T)^\top$. We also let $\Psi=\text{plim}D^{-1}\mathbb{X}^{\top}\mathbb{X}D^{-1}$ and $\Phi=\text{plim}D^{-1}\mathbb{X}^{\top}UU^{\top}\mathbb{X}D^{-1}$, where $D=\text{diag}\big({I^*_1}^{1/2}\mathbb{I}_p,...,{I^*_{m^*+1}}^{1/2}\mathbb{I}_p \big)$ with $\mathbb{I}_p$ denoting the $p\times p$ identity matrix, $U=(u_1,...,u_T)^\top$, and $\mathbb{X}$ is defined analogously to $\widehat{\mathbb{X}}$ but using the true rather than the estimated breakpoints.
\noindentProof: See Appendix D.
We note that Theorem (ref) is a direct consequence of Theorem 3.6 in qian_su_2016 on the asymptotics of their post-LASSO estimator of regression coefficients. We also note that the asymptotic variance matrix can be estimated by replacing $\Psi$ and $\Phi$ with estimates obtained using $\widehat{D},\widehat{\mathbb{X}}$ and $\widehat{U}=Y-\widehat{\mathbb{X}}\widehat{\boldsymbol \alpha}_{m^*}$.
In this section, we investigate the finite sample properties of our proposed approach, focusing on the number of identified breaks.
In order to compare the new estimator to the GFL approach, we follow qian_su_2016 and use the same data generating process as in (ref) with the following cases of interest: the case of no breaks, the case of one break, and the case of many breaks. We compare MIO and GFL methods with classical approaches used by bai2003computation, namely, BIC and LWZ information criteria liu/etal:97, and the sequential method SEQ of bai1998estimating.
BIC and LWZ are natural comparisons as they choose the number of breaks using information criteria. Specifically, they minimize $BIC(m)=\ln{(SSR(\widehat{T}_1,\dots,\widehat{T}_m)/T)}+p^{\star}\ln(T)/T$ and $LWZ(m)=\ln{(SSR(\widehat{T}_1,\dots,\widehat{T}_m)/(T-p^{\star}))}+(p^{\star}/T)0.299(\ln(T))^{2.1}$ over $m$, where $p^{\star}=(m+1)p+m$ and $SSR(\widehat{T}_1,\dots,\widehat{T}_m)$ is the sum of squared residuals from estimating (ref) and the break dates $\widehat{T}_1,\dots,\widehat{T}_m$ are obtained by minimizing $SSR(T_1,\dots,T_m)$ over all possible combinations of $T_1,\dots,T_m$. The commonly used sequential method SEQ is based on first testing for the presence of breaks using the UDmax test of bai1998estimating, and then sequentially testing the null hypothesis of $l$ breaks against the alternative of $l+1$ breaks from $l=1$.
We generate 500 replications for each problem instance described below. For each problem instance, we solve optimization problem ((ref)) in the MIO formulation ((ref))-((ref)). We use a widely available and highly powerful Gurobi 9.5 Solver on one core of AWS EC2 r5.16xlarge (64 cores, 512 GB RAM), with the time limit for each problem set to 450 seconds.
We solve the MIO problem $N$ times for different values of $\lambda \in \{\lambda_1, \lambda_2, ..., \lambda_N\}$ and find the final estimator via the information criterion, as in qian_su_2016. More specifically, we solve $\min_{\lambda_j \in \{\lambda_1, \lambda_2, ..., \lambda_N\}} \log\frac{1}{T}\sum_{t=1}^{T}(y_{t} - \widehat{\boldsymbol \beta}_{\lambda_j,t}^{\top}x_t)^2 + \frac{1}{\sqrt{T}} p (\hat{m}_j+1)$, where $\widehat{\boldsymbol \beta}_{\lambda_j,t}$ is the solution corresponding to $\lambda_j$ and $\hat{m}_j$ is the number of breaks in this solution.\footnote{Python codes and data for simulations and applications are available on the corresponding author's web page and Github.} In practice, we ensure that our grid $\{\lambda_1, \lambda_2, ..., \lambda_N\}$ is sufficiently dense to account for {every} possible number of estimated breakpoints between zero and a computationally reasonable upper bound $B$; for example, $B=25$. We recommend the following approach: (a) if the information criterion is minimized at a $\lambda$ for which $\hat{m}_{\lambda}<B$, then we stop and use $\widehat{\boldsymbol \beta}_{\lambda,t}$ as our final estimator; (b) if the information criterion is minimized at a $\lambda$ with $\hat{m}_{\lambda}=B$, then we increase $B$ to $1.2B$ and re-evaluate the criterion. We continue this process until the minimum of the criterion corresponds to $\hat{m}_{\lambda}<B$.
We report results for the case of one break and many breaks. The case with zero breaks, which confirms our method is competitive in terms of correct detections, is reported in Appendix E. Computational details including boxplots of optimality gaps and run-times are available in an online Supplement.
In this section we use the following DGP variations:
where
and the true parameter values for $\sigma_u$ are the same as in the last subsection.
Table (ref) reports the percentage of replications that detected the correct number of breaks (one) and the accuracy of the break date detection for MIO, GFL, BIC, LWZ and SEQ. Column pce contains the percentage of correct detections, column hd/T reports the Hausdorff distances, divided by $T$, between estimated break date and true break date (conditional on the correct estimation of the number of breaks) and serves as a measure of the overall accuracy of break date estimation.
It can be seen from Table (ref) that MIO and GFL produce comparable results. MIO tends to outperform GFL in terms of $pce$ in larger samples. In fact, when $n=500$, MIO outperforms all the competitors in terms of $pce$, except for DGP-4 in the low signal-to-noise regime ($\sigma_u=1.5$), in which all estimators behave extremely poorly regardless of the sample size. Additional results, not reported for brevity, suggest that MIO tends to overestimate the number of breaks in small samples while GFL tends to underestimate the number of breaks in small samples.
In terms of accuracy of the break date detection, MIO outperforms GFL in all sample sizes but only in high signal-to-noise environments ($\sigma_u=0.5$). In the mid to low signal-to-noise regimes the performance of the five estimators is mixed but MIO shows a superior performance in larger samples.
Overall, the percentage of correct detections decreases for all methods when the noise variance, $\sigma_u$, increases. Both correct detection and accuracy improve as the sample size $T$ increases.
We follow the simulation design of qian_su_2016 and consider
where $x_t \sim i.i.d. N(0,1)$, $u_t \sim i.i.d. N(0, \sigma^2_u)$,
For the first design (DGPn-1), we fix the length of the regime $\Delta=30$ and allow a different number of regimes $R\in\{6,10,20\}$. For the second design (DGPn-2), we fix the number of regimes $R=10$ and allow different regime lengths by varying $T\in\{150,300,600\}$.
Table (ref) reports the percentage of correct detections $pce$ and scaled Hausdorff distance $hd/T$ for each number of breaks $R$ and each sample size $T$. We report $pce$ only for MIO and GFL because the classical methods of bai1998estimating,bai2003computation do not permit 9 and 19 breaks for the selected sample sizes. We can see from the table that MIO uniformly dominates GFL in all cases in terms of $pce$, and the better performance is particularly noticeable in cases of higher error variance ($\sigma_u=0.5$).\footnote{Our results for GFL are different from (and somewhat worse than) those reported by qian_su_2016. We used the regsc package provided by Junhui Qian on his webpage and the same options for all simulations with lambda=NULL and method=“ic” as the defaults. We note that our MIO results dominate both our version of GFL and the original version.} Also, the accuracy of MIO is higher for most cases in terms of the Hausdorff distance.
In this subsection, we consider the U.S. real interest rate time series from 1961Q1 to 1986Q3 used by garcia1996analysis and bai2003computation. Of interest is a simple level shift model which can be written as follows
where $y_t$ is the real interest rate obtained using the U.S. 90-day Treasury bill rate and a quarterly inflation rate series constructed from the U.S. CPI.
Table (ref) reports the estimation results using the five methods. The breakpoints detected by MIP are depicted in Figure (ref). The MIP and GFL methods detect 4 breaks while the classical methods BIC, LWZ and SEQ detect 2, 0 and 3 breaks, respectively. The most common break date is 1972Q4 which can be tied to the oil crisis, while the breaks in the 1980s can be associated with Paul Volcker's attempt to lower the inflation in the US and with the start of what became known as the Great Moderation, that is, a decrease in the variance of the error term. For example, walsh1987testing noted that in early 1980, real interest rates may have reacted to a change in the Federal Reserve operating procedures implemented in late 1979.
We can see that the level breaks detected by MIP and GFL in the 1980s are very close to one another. This questions the assumptions of BIC, LWZ and SEQ about the minimal length of the regime, which may affect these methods' performance.
Our finding of the additional break in 1983Q1 provides support to bauwens/etal:15 who argue that a distinct regime starting in 1983Q1 is characterized by the lowest variability, the highest persistence and an average growth rate between the other regimes, and thus best defines the start of the Great Moderation. This finding also contributes to the literature on what caused the high real interest rates in the U.S. in the early 1980s. Specifically, it puts more weight on the Federal Reserve's attempts in late 1982 to de-emphasize monetary aggregates, described by walsh1987testing. As argued by huizinga1986monetary, this suggests that it is the monetary policy, not fiscal deficits, that bears the responsibility for high U.S. real interest rates in the early 1980s.
In this subsection, we consider a business statistics application from inventory modeling. A classic model in this literature is \[z_t = \mu +(1-\alpha) z_{t-1} + \epsilon_t, \] where $\epsilon_t$ is white noise, $z_t$ is a disequilibrium variable, e.g., inventory-to-shipment ratio, and $\alpha\in[0,1]$ is the speed of adjustment caballero/engel:93, jorda1999random, jorda/marcellino:04. A problem with explaining the inventory dynamic is that the empirical values of the adjustment speed are much lower than predicted by economic theory. For example, estimates of $\alpha$ referenced by jorda/marcellino:04 and jorda1999random are between 5% and 29% for monthly data suggesting that it takes over a year to close 95% of a given disequilibrium.
One solution proposed by jorda1999random is to recognise that data aggregation happens at irregular intervals -- different from the reporting frequency -- and affects $\alpha$. For example, jorda/marcellino:04 estimate a Markov switching model with three regimes representing what they call zero, one and two original-time adjustments, that is, three types of time-varying aggregation frequency. The values of $1-\alpha$ they find are 0.99, 0.72 and 0.58. Because of the limitation on the minimum number of observations in each regime and the small sample, their number of regimes is fixed at three.
We re-estimate the model without that restriction. We use the data on total manufacturing from the Census Bureau’s monthly survey “Manufacturers’ Shipments, Inventories and Orders”. This is the same data source as used by jorda/marcellino:04 and jorda1999random but a larger product category (they used glass containers which are no longer available separately), nominal volumes and a larger sample ranging from January 1992 to May 2023. Figure (ref) shows the inventory-to-shipment ratio for the entire sample of 380 observations. Similar to jorda/marcellino:04, the data is not seasonally adjusted, which is preferable when investigating inventory dynamics.
When we apply our method to this data, we find four regimes, which are delineated in Figure (ref) by the dashed lines and correspond to January 1997, October 2008 and November 2019. The estimated values of $(1-\alpha)$ are 0.142, 0.135, 0.094, and 0.041, suggesting the adjustment speed $\alpha$ of 85.8-95.9%, which is much faster than in previous studies. A simple calculation used by jorda1999random suggests that 95% of a disequilibrium is replenished within $T=30 \log(0.05)/\log(0.142) =46$ days, i.e., within a month and a half -- a more realistic estimate.
For comparison, the GFL estimator produced zero breaks; BIC and SEQ resulted in four breaks but the estimates of $(1-\alpha)$ include negative values, which are infeasible and, thus, not reported. LWZ produced three breaks in July 1997, October 2008 and May 2019, and the estimates of $(1-\alpha)$ are 0.195, 0.140, 0.133 and 0.036. Three of these estimates are very close to ours while the highest value of 0.195 means that the longest period required to close 95% of a disequilibrium is about two months, i.e., only two weeks longer than our estimate.
jorda/marcellino:04 hypothesize that these regimes reflect time aggregation, which in effect scales the original sampling frequency. If this is the case, one can make an argument similar to jorda/marcellino:04 that the estimates of $\alpha$ corresponding to the different aggregation scales should relate to each other. A smaller estimate of $1-\alpha$ is expected to be close to a power of a larger estimate. It turns out that our largest and smallest estimated values of $1-\alpha$ (i.e., 0.142 and 0.041) agree with this prediction: $0.142^2 = 0.02$, which is not far from 0.041, or, equivalently, the implied $\alpha$ derived from $(1-\alpha)^2=0.041$ is 0.8, which is not far from 0.858. We note that a similar relationship exists in the LWZ estimates: the implied $\alpha$ derived from $(1 -\alpha)^2 = 0.036$ is 0.81, which coincides with $1-0.195=0.805$. The intermediate values of $1-\alpha$ found by our method and by LWZ do not seem to obey this relationship, suggesting that these regimes may have a different nature. Nonetheless, the above findings can serve as evidence in favor of the hypothesis of jorda/marcellino:04 for much higher speeds of adjustment than previously thought.
We propose a new way of handling change-points in econometrics based on computational advances in mixed integer optimization and we work out statistical properties for the estimator of the number of breaks, break locations and the regression coefficient in one step. The approach shows remarkable adaptivity and versatility in that it has similar or better asymptotic properties than the LASSO-based alternatives under somewhat weaker assumptions. In simulations, the proposed method shows improved performance especially under many breaks scenarios. In empirical applications, it provides additional insights, offering a new and robust way to obtain evidence on the number and location of breaks, and the corresponding regime characteristics without the restrictive assumptions on the regime duration. We demonstrate this for two classic examples in economics.