The problem of detecting structural breaks in time series is well-recognized in the literature. There are two base variants: online and offline, the use of which is determined by the practical setting in which we need to execute this task. There were many works devoted to the online variant such as these by Tartakovsky (2019) or Huang et al. (2020). However, the online processing is outside the scope of our interest. We are only focusing on the offline processing mode.

Offline structural break detection has its own subvariants. The most classical division of available methods is into two groups. The first group aims at finding at most one changepoint (AMOC) detection, where the aim is to state if a time series contains one structural break and, subsequently, where it is located. The popular methods include cumulative sums (CUSUM) (Yang and Zhang, 2022), (SCUSUM) (Shi et al., 2022b) and likelihood ratio tests (LRT) (Bai et al., 2024). Nonetheless, in this article we will focus on multiple structural breaks detection problem.

We can also divide the existing structural break detection methods according to how the structural break is defined and found. The most straightforward approach is a canonical segmentation problem, where the only parameter analysed during breakpoint detection is time series mean. Numerous works have been dedicated to this approach, and they have been thoroughly discussed by Cho and Kirch (2024). The more advanced approaches consider changes in time series variance or model data with the use of one of the common time series modelling methods. When it comes to the latter aspect, we shall underline that different models may be used for this purpose, for example, the autoregressive model (AR) (Behrendt, 2021) or the GARCH model (Cho and Korkas, 2022). Many papers are dedicated to simple and very narrow problem modelling. For instance, Romano et al. (2022) use AR(1) model.

This article is devoted to multiple structural break detection in an offline manner, so let us describe analogous leading and popular methods used for the task. As Shi et al. (2022b) say, we can distinguish two categories of methods for solving the task. First are methods using recursive partitioning of a time series into smaller segments or other approaches utilizing AMOC methods. Second are methods analysing the entire time series.

A distinct category can be assigned to the methods which sequentially analyse time series and are able to increasingly add new structural breaks, such as the method proposed by Bai and Perron (1998), which is still recognized as the baseline method in economics (for example, Lee et al., 2022). This approach, commonly known as the Bai-Perron method will be mentioned later in this paper.

The first mentioned family of approaches (recursive one) generally uses AMOC methods to detect one structural breaks and later run the same AMOC method for chunked data. The base method in the area is binary segmentation (BS) (Scott and Knott, 1974). It was later extended by Fryzlewicz (2014), who proposed wild binary segmentation (WBS). There are still works in this direction and new modified variants of binary segmentation were being proposed, such as WBS2 (Fryzlewicz, 2020) or seeded binary segmentation (SBS) by Kovács et al. (2023). Another interesting approach of this kind is Classification Score Profile (ClaSP) presented by Ermshaus et al. (2023).

Many different approaches were proposed that can be described as general methods that jointly analyse entire time series. Some are using Bayesian learning or other probability-theory-based techniques for changepoint detection. In this line of research, the algorithms typically employ recursion to partition the data into segments. We envision that the changepoints of interest separate these segments. Bayesian approaches use different methods to solve the problem. Another technique was delivered by Li et al. (2019a) using a linear regression model, AR in particular. Yet other studies use the Markov chain Monte Carlo sampling, specifically the Metropolis-Hastings algorithm, to search for structural breaks. Shaochuan (2020) studied the use of hidden Markov models with continuous-time FFBS and continuous-time Viterbi algorithm. Casini and Perron (2024) proposed a Laplace-based (Quasi-Bayes) procedure utilizing the Monte Carlo method.

Another notable way to solve the problem of locating multiple structural breaks in a time series is Cumulative Sums of Squares (CUSUM), studied by Inclan and Tiao (1994). Among more recent studies considering the CUSUM approach, we can mention the works of Madrid Padilla et al. (2022), Kim et al. (2022) or Borzykh and Yazykov (2020). The CUSUM approach passes through the time series to search for points where the sum of squares of previous and present values changes.

Other approach to the problem was moving sum (MOSUM) procedure. Their method is essential in canonical segmentation problem. It was recently used, for instance, by Cho and Kirch (2022) or Meier et al. (2021). MOSUM was utilized also for more advanced segmentation problems such as in the work of Kirch and Reckruehm (2024).

It is vital to notice that many of known approaches use constraint on the minimum distance between two consecutive structural changes. The approach is very popular and was used, for example, by Davis et al. (2016), Doerr et al. (2017) or Yan et al. (2021).

As a side note, we have observed a growing number of papers elaborating on the benefits of including a structural break detection step in time series forecasting procedures. Among researchers focusing on this aspect, we find Smith (2023), Safikhani et al. (2022), and Altansukh and Osborn (2022). Another substantive stream of studies emerged when the existing changepoint detection methods were adapted to temporal panel data, as indicated by Bardwell et al. (2019).

The last group of methods for solving the task we want to discuss are evolutionary algorithms. In particular, it is worth mentioning the works of Davis et al. (2016, 2006, 2008), Davis and Yau (2013), who conducted several studies on the application of the genetic algorithm in the discussed context. They have referred to their new method as Auto-PARM. The chromosome in this algorithm contained information about whether a structural change was present at a given location. The method uses relatively simple mutation, various crossover variants, and proportional selection. In addition, the island model of the genetic algorithm was employed to improve the results.

The genetic algorithm was also considered in the works of Doerr et al. (2017). As already mentioned, their algorithm was constructed for long time series up to one million data points. The results of this study indicated that using both single-point and uniform crossover simultaneously improves the quality of the results. Furthermore, the effectiveness of enhancing the algorithm by using multi-start and champions league strategies was noted. A similar outcome was documented by Suárez-Sierra et al. (2023). On the other hand, Lim et al. (2020) presented a memetic algorithm for multivariate time series segmentation dedicated to long time series.

Éltetö et al. (2012) proposed using the highly acclaimed CMA-ES (Covariance Matrix Adaptation Evolution Strategy) algorithm. They divided their method into two steps. The first step aims at finding the number of structural changes. The second step optimizes the locations of each structural change. They did this by analysing solutions for different, algorithmically chosen numbers of breakpoints. CMA-ES operates in a continuous space, while the change point detection problem is embedded in a discrete space. Therefore, it was necessary to move the problem to another space, which was done by processing the structural change locations and rounding them for solution evaluation.

Many of the methods used for structural break detection require an objective function that allows comparing particular structural change positions. Choosing an appropriate objective function is not easy because of the need to simultaneously fit the data model and minimize the number of model parameters, as state (Shi et al., 2022b). Different objective functions were used in the literature. Some were very simple and failed to reflect the nuances of model fitting. This was due to various reasons. In some cases, the computational complexity was decisive, as it was in the approach of Doerr et al. (2017).

Penalty function methods have been analysed as well. For example, Hall et al. (2013) used a penalty function in which three parameters of a time series segment model corresponded to one structural change. Another variant of the penalty function was considered by Doerr et al. (2017), whose penalty function was the sum of the penalty for the number of detected breakpoints and the penalty for the fitting error calculated separately for each segment of the time series as the area of a rectangle whose sides were: the beginning and the end of the studied segment, and the minimum and maximum value of the series on the given segment.

For more complex approaches, a fairly obvious choice was to check known information criterion models. Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), or modified Bayesian Information Criterion (mBIC) can be employed (Zhang and Siegmund, 2007). An intensive comparison of these criteria and MDL (see next paragraph) was delivered by Shi et al. (2022b).

In recent years, the principle of Minimum Description Length (MDL) formulated by Rissanen (1978) has gained substantial popularity in this field. It was applied to the problem of structural change detection by Davis et al. (2006). This solution was originally adapted for use together with the time series autoregressive (AR) model. Later works made it possible to transfer this method to other time series models: Generalized AutoRegressive Conditional Heteroskedasticity (GARCH) (Davis et al., 2008), periodic AR (Lu et al., 2010) or AutoRegressive Moving Average (ARMA) (Davis and Yau, 2013). As mentioned, Li et al. (2019b) formulated a modification of the MDL method known as the Bayesian Minimum Description Length rule. The cost function based on the MDL principle was later used in many different works, such as by Ditzen et al. (2021) or Lan et al. (2022). A modification of the MDL rule was also shown in Bayesian Minimum Description Length (BMDL) method proposed by Li et al. (2019a).

Time series can be simply described as a sequence of values observed in time. There are many different ways to model time series. One of the very popular is the use of models from the AutoRegressive Integrated Moving Average (ARIMA) family. The simplest model belonging to this family is the autoregressive model (AR), which can be written as (1) y t = c + β 1 y t − 1 + β 2 y t − 2 + ⋯ + β p y t − p + ε t , where: •

c, β 1 , … , β p – model parameters,

The model can be briefly expressed as ARIMA ( p , d , q ), where p is the order of the autoregressive part, d is the degree of differencing, and q is the order of the moving average part.

A time series structural break is a place in the time series where data changes its characteristics. It can be a change in the time series mean, variance, or a point where a trend has emerged or disappeared. The most often-mentioned definition says that it is a place in the time series where the time series model or one of its parameters changes (for example, Cheng et al., 2022). Obviously, time series can have several structural breaks or not have them at all.

In this work, we use metaheuristic optimization-based algorithms to solve the problem of time series structural break detection. We selected Ant Colony Optimization (ACO), Particle Swarm Optimization (PSO), and Genetic Algorithm (GA). The algorithms mentioned are all metaheuristics, which means that they do not guarantee that an optimal solution will be found. However, they can find a good-enough solution in a relatively short time. This property makes them indispensable in many engineering applications, see Abdel-Basset et al. (2018).

ACO, PSO, and GA are based on an analogous idea that the algorithm operates on a collection of candidate solutions called population. The procedure is iterative. In each iteration, the population undergoes specific modifications, and in this manner, candidate solutions are constructed or modified. An essential part of these algorithms is a cost function that evaluates how good a particular candidate solution is. Since these are heuristic optimization approaches, termination criteria concern either a count of performed iterations or a count of iterations without any improvement.

ACO is a method proposed and developed by Dorigo and Stützle (2019). It is patterned on the behaviour of ants. The population of the algorithm is interpreted as a swarm of ants that create solutions by walking on a graph. All solutions found in a given iteration are evaluated, and ants put pheromones on the edges of a graph used in these solutions according to their fitness.

Ants construct solutions in a probabilistic, step-wise manner. If we have an iteration t and an ant, which built a part of a solution and is in a vertex i, the probability of the ant’s move from the ith to the jth vertex is given as (4) p i j k ( t ) = ( τ i j ( t ) ) α · η i j β ∑ a ∈ N ( τ i a ( t ) ) α · η i a β , j ∈ N , 0 , j ∉ N , where: •

τ i j ( t ) – pheromone count between vertices i and j at the moment t,

The general scheme of ACO is presented in Algorithm 1.

PSO is patterned on the behaviour of birds’ flocks or schools of fish. In contrast to ACO, the algorithm operates in a continuous space. The population is interpreted as a collection of particles. Each particle has its own location, memory, and velocity. The location of a particle i is denoted as x i , the velocity is denoted as v i , the best solution found by the particle is denoted as s i , and the best global solution is denoted as s b . During each iteration the velocity and location of each particle are modified according to the formula: (5) v i ( k ) = ω v i ( k − 1 ) + c 1 r 1 [ s i ( k − 1 ) − x i ( k − 1 ) ] + c 2 r 2 [ s b ( k − 1 ) − x i ( k − 1 ) ] , where:

The location of a particle is updated according to: (6) x i ( k ) = x i ( k − 1 ) + v i ( k ) .

The general flow of the PSO is presented in Algorithm 2.

GA is a long-known evolutionary algorithm that uses a population of solutions. Each solution is represented by a chromosome, which consists of genes. The population is modified by three fundamental operations: selection, crossover, and mutation, which are performed during each iteration. These operations are defined very generally. Thus, so many variants of the GA have arisen. The first operation is selection, which aims to pick the best individuals from the population. The two most common selection mechanisms are proportional selection and tournament selection. A mutation usually makes a minor change to the solution. For example, in the popular bit flip mutation for binary encoded chromosomes, one or more bits are randomly selected and inverted. There are other types of mutations, such as swap mutation, truncation mutation, or inversion mutation. Crossover initially results in several parent solutions that combine to form one or more offspring. Popular crossover methods are uniform crossover, one-point crossover, or multi-point crossover. In many situations, applying standard mutation or crossover methods is impossible, and the need to modify them or create custom, problem-specific crossover or mutation operators arises.

The general scheme of the GA is presented in Algorithm 3.

The island model is one of the approaches to perform multi-population optimization. In the model, there are several populations located on several islands. After every t iterations, we perform a migration that allows replacing the worst solutions from one population with the best solutions from another population. The count of iterations between migrations, population size, and the count of migrated elements are parameters of the model.

The island model is a simple, yet effective method to get good results fast. It is often used with the genetic algorithm but can also be applied to other evolutionary algorithms.

Selecting a cost function to measure the quality of solutions is a critical problem in metaheuristic optimization. Constructing a suitable cost function for the structural break detection problem is difficult because we need to minimize the model errors and the number of structural breaks and balance the two.

The literature offers various solutions to this problem; more information is available from Farsi et al. (2020). One of the most acclaimed ones was the use of the minimum description length (MDL) rule. This rule states that a model is only as good as long as it needs little space to save its data. We can write it as (7) C L ( y ) = C L ( M ) + C L ( e | M ) , where C L is a code length, y is the analysed time series, M is the model describing it, and e | M is the error of fitting M to y. C L ( M ) is simple to calculate compared to C L ( e | M ). Rissanen (1978) proposed to calculate it using the log-likelihood function: (8) C L ( e | M ) = − log L ( M ) .

Davis et al. (2016) introduced an adaptation of the MDL technique to the problem of structural breaks detection for the AR model: (9) MDL ( m , p 0 , ( τ 1 , p 1 ) , … , ( τ m , p m ) ) = log 2 m + m log 2 n + ∑ i = 0 m log 2 p i + ∑ i = 0 m p i + 2 2 log 2 n i + ∑ i = 0 m n i + n i log 2 ( 2 π σ i 2 ) 2 , where: •

m – structural change count,

Let us define the MDL function for the ARIMA model based on Eqs. (7) and (8) as (10) MDL ( m , ( τ 1 , p 1 , q 1 ) , … , ( τ m , p m , q m ) ) = log 2 m + m log 2 n + ∑ i = 0 m log 2 ( p i + d i + q i ) + ∑ i = 0 m p i + d i + q i + 2 2 log 2 n i − max i log L ( M i ) , where log L ( M i ) is the log-likelihood function.

The defined function in Eq. (10) relies on the addition of four components. The first component ( log 2 m) represents the number of structural breaks. It is a natural number. The logarithm is used to represent this value in terms of bits, according to the idea of the original MDL function. The second component is m log 2 n. It represents the space needed to save the location of each structural break. Each location is kept as a natural number that is not larger than the series length n. The third component represents the ARIMA specification ( p , d , q ) for each segment in the time series. The final component keeps the values of ARIMA model parameters for each segment. These values are computed based on n i elements.

In this subsection, we address new cost functions. Each of them was defined in two variants – for the AR and ARIMA model.

We proposed two extensions to the cost functions that can increase their capabilities.

The first extension is a weight for enforcing minimal allowed interval length. A common pitfall of structural break detection algorithms is that they find too short, insignificant time series segments. So far, the problem has mostly been solved by inserting a strict ban on generating time series segments shorter than a certain threshold. Such a tough, binary solution is not clear and justified. Thus, we propose a new attempt to solve this problem. We suggest adding to cost functions a penalty term defined as follows: (21) P F S ( s i ) = ( b − n i ) 2 , n i < b , 0 , n i ⩾ b , where: •

P F S ( s i ) – function value for a segment s i ,

The second extension was dedicated to the methods that do not account for trends. We propose to solve this deficiency by adding a weight to a cost function to estimate the deviation from a linear trend. We can define it as (22) M S E ( s i ) = ∑ j = 0 n i − 1 ( s i j − a j − b ) , where

Both extensions can be applied with the use of additional parameters, giving them suitable impact.

While designing the new solutions, we paid a lot of attention to the issue of computational complexity. It was critical since the problem at hand involved intensive computations. The algorithms we used needed information about cost function value linked with time series segments. The information could be calculated dynamically during the computation of a structural break or eager (earlier) before the computation. The second strategy seems more straightforward as we divide our task into two sub-tasks: computation of a cost function for all possible segments and finding structural breaks. It is viable only for short time series because executing it for long time series would induce many redundant calculations.

Let us address in more detail the computational complexity of finding structural breaks. We will mark population size as P, iteration count as N, length of time series as n, and time needed to get cost function for the segment as C. If we use precalculated values for it, C is constant. In PSO, two procedures: modification of the solution and its evaluation, are performed in each iteration. The time complexity of modification is proportional to the maximum count of structural breaks in the solution, which is a parameter of the algorithm, and we will denote it as K. So, the time complexity of evaluation is proportional to the product of K and C. Thus, the time complexity of PSO can be expressed as: (25) t PSO = O ( P · N · K · C ) .

In each iteration, GA performs selection, mutation, crossover, and evaluation. The time complexity of two-point crossover and mutation is O ( n ) and the time complexity of tournament selection is O ( T ), where T is a tournament size. The time series evaluation cost is K · C. However, K is not constant, but it is a number not larger than n. Thus, the time complexity of GA is given as: (26) t G A = O ( P · N · ( T + n · C ) ) .

ACO performs three operations: construction of a new solution, evaluation, and updating pheromone trace in each iteration. During solution construction, an ant calculates probabilities and needs knowledge of cost function at analysed time series segments, which is done in O ( n · C ) time. The maximum count of detected structural breaks is O ( n ). Thus, the construction of a single solution is calculated in O ( n 2 · C ) time. The solution evaluation is in O ( n · C ) time. The update of the pheromone trace is made for all ants together and has three important parts: 1.

Sorting solutions calculated by ants in time O ( P 2 ).

In this section, we present the experiments’ outcomes to validate the quality of the proposed novel methods. The scope of the study covered both artificial and real-world time series representing different situations. Most of the analysed time series were relatively short, since conventional structural break detection tasks are executed for such cases. We analysed eight synthesized time series, which included instances with different numbers of structural breaks, appearing and disappearing trends, and changing time series mean and variance. There was also an example of a white noise time series with added variation, in which the algorithms should not detect any structural breaks. We have defined “test cases” of synthetic time series. We assumed some underlying time series properties for each test case and then generated each test case multiple times, maintaining the same underlying properties but with some random noise. This entailed that the experiments were repeated multiple times for each test case. Specific differences between particular time series instances within a single “case” arise only because we draw values, but the parameters of distributions we use stay the same.

On top of ten synthetic series, we used two real-world time series representing the power demand of a fridge freezer in a kitchen. They come from the Electrical Load Measurement dataset by Murray et al. (2015).1¹

The dataset is available at https://www.timeseriesclassification.com/description.php?Dataset=FreezerRegularTrain

Let us list and briefly discuss synthetic test cases:

To eliminate the possible bias arising due to the randomness of the empirical methodology of this study, we have been performing multiple repetitions of the experiments. For each synthesized test case, we generated 10 time series instances with different random number generator seeds. Experiments were executed for each generated time series. As a result, we performed 10 repetitions of each experiment for each test case. Our motivation was to ensure a reasonable variability of the data by maintaining the same data generation process that relies on drawing from distributions. To keep the same level of objectivity in result evaluation, in the case of the real-world time series, we were also repeating the experiments 10 times.

The outcomes of the experiments that were conducted involving different algorithms were compared with each other and with the outcomes obtained by other researchers. We assume three state-of-the-art methods, which we use for such comparisons. The first is the method proposed by Bai and Perron (1998), Bai and Perron (2003), and the second is ClaSP by Ermshaus et al. (2023). The third is the method by Davis et al. (2016), which can be seen as the immediate predecessor of the approach introduced in this paper. In the first two cases, we run the algorithm to produce the results, but in the third case, we compare our results to those presented by Davis et al. (2006) using analogous time series examples.

We propose two metrics for the evaluation of the outcomes of structural break detection procedures: •

break point difference (BPD) – absolute value of the difference between the found and correct count of structural breaks,

At first, we compared results obtained with the MDL, SPF1, SPF2, EMDL, DCF, and PDCF cost functions. We used both the AR and ARIMA models and both versions of the cost function: with and without a weight for deviation from a linear trend. The experiments concerned seven synthetic time series and two real-world time series. The only omitted case is Fig. 1(h), which will be discussed in detail in Section 7.

Experiments were performed using PSO with parameters C 11 : ω = − 0.1832, c 1 = 0.5287, c 2 = 3.1913, population size 47. These parameters come from Hvass Pedersen (2010). We have tested 188 variants of cost functions: 164 for the AR model and 22 for the ARIMA model. In principle, we have employed a grid search of procedure parameters for SPF1, SPF2, EMDL, DCF, and PDCF functions with the AR model. In the case of ARIMA, we have performed fewer tests due to a more considerable computational complexity – we tested a subset of configurations close to the configurations, for which results for the corresponding functions with the AR model were good.

A detailed discussion on optimization algorithms hyperparameters and their fitting is provided in Section 6.

Our results were evaluated with the use of BPD and SC metrics and conclusions attained with use of both were very similar. The best results (assessed with the use of BPD and SC metrics) were obtained consecutively for configurations: EMDL AR ( μ = 3, ν = 1), EMDL AR ( μ = 3, ν = 1), MDL ARIMA , all without involvement of the MSE modification. There were also other cost functions which led to noteworthy results. The most interesting configurations are presented in Table 1 and Table 2.

Let us discuss the behaviour of various cost functions in more detail. For the simplest cases, analysed functions perform very well. For test case 1(a) most of the cost functions in almost each run correctly detect places where the average value of time series changes. The exception was SPF AR ( μ = 0.1, ν = 1) with no MSE modification, which did not find any structural break.

Most of analysed methods correctly find the place where the variance changes (test case 1(b)). Problems have only SPF and SPF2 methods. They generally did not detect a break in such a case. The SPF AR ( μ = 10, ν = 10) with MSE is an exception, because it finds structural breaks but locates them in incorrect places.

Most of presented cost functions detect change between two different AR models very well (case 1(c)). The exception are methods involving MSE modification. They detect usually two structural breaks and only small part of them is close to correct structural break location.

Most of presented cost functions correctly solved the task of processing white noise time series (case 1(d)). The expected outcome in this scenario is not to return any structural breaks. The only cost functions which almost always failed the test are the PDCF and DCF cost functions.

Test case 1(e) contained series composed of two segments: white noise and a segment with a trend. In this scenario, the classical MDL and the elastic MDL rule produced excellent outcomes. PDCF and DCF cost functions were less precise, but they detected one structural break near the expected location (slightly too early). Methods based on penalty cost function in most cases returned more than one structural break. Usually one of them was in the expected location (with the exception of configuration SPF AR ( μ = 0.1, ν = 1) without MSE where no structural breaks were found.

Result obtained for the test cases composed of three segments of the same length (case 1(f)) were generally good. Only SPF and SPF2 cost functions had problems. The configurations did not detect any structural break or detected (often more than one) in incorrect places.

The next test case concentrated on a rather theoretical example, where time series was the shape of a parabola (see Fig. 1(g)). The parabola case is very difficult as there are many viable locations of structural breaks. In this test case, the introduced BPD and SC scores are not very informative and a valid evaluation must be manual. The experiments have shown, that the methods including the MSE modification and SPF2 AR ( μ = 0.1, ν = 10) detected too many structural breaks and were not very practical. Conversely, cost functions employing the ARIMA model had a strong tendency not to detect any structural change. Results generated by SPF AR ( μ = 0.1, ν = 1) also were very impractical (one structural break detected closely to start or end of time series). Other cost functions usually allowed for detecting a few structural changes near the 1 4 th and the 3 4 th of the time series length. The best results were obtained with the EMDL AR ( μ = 3, ν = 1) model.

The last two examples were the real-world time series describing fridge freezer power demand. We analysed two examples. First, where the freezer was turned on and turned off after some time (case 2(a)). Second, where the freezer was turned on and was working until the end of the time series (case 2(b)). The best results were obtained for chosen configurations of SPF and SPF2 methods. For case 2(a) they allowed for a correct detection of both structural breaks and for case 2(b) the same methods returned one structural break in the correct location. It is worth to say that cost function variants using MSE modification give more reliable results for case 2(a). All other cost functions listed in the tables did not find any structural break or found one near to the end of the time series. The latter behaviour is a result of a phenomenon which will be discussed in the next subsection.

We have observed that the proposed cost functions enhance time series processing capabilities for structural break detection algorithms. Noteworthy, for all of the proposed cost functions, we can find a set of parameters that works very well for various test cases. Furthermore, we will demonstrate in Section 7 that the proposed cost functions are generally better performing than the state-of-the-art MDL AR cost function. Promising results have been achieved with the use of the SPF and SPF2 cost functions. Their superiority was the clearest when we processed real-world time series. Their only weakness is that they may tend to output incorrect structural break count – too many or too few, more often too many, especially when the methods are used with MSE modification.

The empirical experiments have shown that the cost functions using the ARIMA model are very good, and in many cases, they are better than the ones using the AR model. Their weak point is computational complexity, which is higher, and this difference grows as the time series length increases. For example, it took only 0.3s for the AR model to calculate MDL cost function partial values for all possible time series segments for the test case with 50 data points and 190.0s for ARIMA. In contrast, for a case with 100 data points, it took 1.4s and 937.7s, respectively. This shows that cost functions using the ARIMA model are helpful for short time series. In practice, it is not a severe limitation for many real-world domains applying break detection methods (such as economics), in which relatively short time series are of interest.

The proposed objective functions punish for deviations from theoretical models constructed on developed segments. Therefore, if these models cannot reflect the underlying data well, there is no reasonable justification to apply this model. If the data contains a clearly visible and changing trend, the recommended choice is an objective function that punishes deviations from theoretical models. The component that punishes for short intervals should we weighted as less important. Punishing for deviations from a theoretical model is a lousy strategy for high-variance time series or if variance is the crucial distinction for subsequent segments. In this case, punishing for short segments and diminishing the punishment for deviations from theoretical models is much better.

We performed tests for 3 configurations of PSO:

All the experiments were performed with maximal limit of 10 detected structural breaks. Tables 3 and 4 present results for these configurations.

According to the two presented metrics, the best results achieved using PSO were for the parameter configuration C 11 . All tested configurations produced excellent results for the first six case studies (1(a)–1(f)).

The most interesting observations were made for the real-world time series. All named configurations did not locate structural breaks in the correct locations and found none or one structural break close to the end of the time series. The detection of such structural breaks (near or precisely at the end of a sequence) is a side effect of the algorithm’s operations. During the optimization procedure, each structural break location is treated as a variable, and it can change its position in time series. A variable at the end of the time series is treated as if the structural break did not exist. This, however, may result in the observed adverse effect. If a time series is challenging, we may end up with incorrectly identified structural breaks near the end. The described behaviour did not occur for configuration C 13 , but it appeared for C 12 , where the population count was smaller (20), and for C 11 , probably due to the negative inertia.

For the GA we ran the experiments on a grid of parameters. We have tested mutation probabilities from the set { 0.0125 , 0.0083 , 0.0033 }, crossover probabilities from the set { 0.2 , 0.5 , 0.8 }, tournament sizes { 2 , 8 }, population size 40 and using small direction modification (0.2) or not.

Both metrics indicated the same best-performing parameter configurations listed below:

The results calculated for listed configurations are in Tables 5 and 6.

The analysis showed that it was generally better not to use the direction modification. The average value of the SC metrics of the results concerning the direction modification was 1.467, while the average for the results without it was 1.187. For various values of mutation probability, crossover probability, and tournament size, mean SC metrics fluctuations were very little, especially compared to concerning mutation modification influence.

Qualitative analysis of the experiments’ results showed that the outcomes were quite similar for all tested GA parameters. There were minor differences in the count of detected structural breaks in a few cases.

The experiments with ACO were conducted for a grid of parameters. We tested different combinations of parameters assuming the following sets of values of interest: the evaporation parameter { 0.5 , 0.9 }, importance of pheromone count { 1 , 1.5 }, importance of fragment rank {2, 5}, and estimated break count { 1 , 2 }, count of elitist ants used to correct pheromone data { 2 , 4 }, importance of count of pheromone updated at one time { 0.05 , 0.2 } and the limit of detected structural breaks { 4 , 10 }. The best results were obtained for the following configurations: •

C 32 : evaporation parameter: 0.9, pheromone count parameter: 1, fragment rank parameter: 2, estimated break count parameter: 2, elitist ants count: 4, updated pheromone count: 0.05, detected structural break limit: 4,

Subsequently, we looked into each parameter and analysed its role in the achieved outcomes. As BPD and SC metrics usually give similar conclusions, but the SC metric provides more information, the discussion below refers to the average value of the SC. The results showed that:

The most significant conclusion refers to the structural break count limit. This limit allows for accelerating calculations and achieving better results for test cases where the desirable structural break count is less than this limit. The experiments showed that more favourable results were obtained with a bigger limit, which indicates that the algorithm can very well generalize and detect fewer breaks than this limit, which is an essential feature for many possible real-world applications of structural break detection methods.

We performed experiments using the Bai-Perron method with the same test cases as in previous tests. Mean BPD and SC for these examples and general mean are presented in Table 9.

We can see that the Bai-Perron method was substantially worse than the presented methods. Qualitative analysis showed that the Bai-Perron method achieved moderately satisfying results for cases 1(c), 1(e), or 1(f). However, even in these cases, the structural breaks were inexact. The Bai-Perron method produced even worse outcomes for the remaining test cases, confirming results for individual test cases. It is worth highlighting that the method had problems with real-world cases where it detected too many structural breaks. The experiments proved that the methods presented in this paper guarantee better results than the baseline Bai-Perron method (please compare with Tables 1 and 2 concerning our results).

We compared our work also to the state-of-the-art ClaSP algorithm by Ermshaus et al. (2023). We run the algorithm on our examples. The obtained results were very poor. They are presented in Table 10. The method found any structural breaks rarely. The only example from set 1(a)–2(b) where structural break was detected was example 1(c), but the precision also was not ideal.

The results covered in Table 10 (ClaSP) may be compared with Table 1 and 2 concerning our results. The advantage of the proposed approach is visible.

Contrary to previous experiments, to process the real-world dataset published by Davis et al., we used the island model (40 islands, 20 migrations, 5 iterations between migrations, small population size 40) to make our settings similar to the settings employed by Davis et al. (2005).

Table 11 presents mean precision of detected structural break locations and associated standard deviations. Particular hyperparameter configurations used with our optimizers are noted with letter C, and their explanation is provided in Section 6.

Davis et al. (2005) delivered a method that returned two structural breaks for every experiment. Our algorithms also found two structural breaks for every performed test.

Subsequently, we compared our results against those calculated using Bai-Perron’s and ClaSP methods. They are also given in Table 11. Please note “Std” values, which show the efficiency of the compared methods (the smaller the better).

It is essential to state that the means and standard deviations presented in Table 11 are calculated only for experiments where two structural breaks were detected. Notably, Davis’s algorithm and our configurations consistently detected two breaks each time we ran them. That was the expected behaviour. However, the Bai-Perron method detected two breaks in 80% of repetitions, while ClaSP detected two breaks only in 40% of repetitions. That information is sufficient to claim that both Bai-Perron and ClaSP methods work worse than the other analysed methods.

Still, to provide a deeper analysis, Table 11 presents the results achieved by all four approaches limited to the cases when two structural breaks were detected (which promotes Bai-Perron and ClaSP, as we skipped the cases when these two returned other breaks numbers).

First and foremost, Table 11 allows comparing our results and those of Davis et al. In both cases, two (out of two) structural breaks were always detected. However, using our method, the second structural break was detected much better than in Davis’s work, while the first was slightly worse. If we take an average of standard deviations for the location of the first and the second structural break, we see that the average for Davis’s work is 0.0045, and the average of the best variant of our method is 0.00212. The latter average (our approach) is two times better than the Davis’s. This shows that the method presented in this article outperforms the one given by Davis et al.

Our approach is very good in comparison with state-of-the-art methods. In our opinion, it is also more versatile and customizable. We have delivered not a single method but a suite of methods. Finally, we may emphasize that our results are much better than those obtained with Bai-Perron’s and ClaSP methods. These two algorithms struggle to find two structural breaks. Even if they find two structural breaks, they locate them inaccurately.