Ensemble binary segmentation for irregularly spaced data with change-points

We propose a new technique for consistent estimation of the number and locations of the change-points in the structure of an irregularly spaced time series. The core of the segmentation procedure is the ensemble binary segmentation method (EBS), a technique in which a large number of multiple change-point detection tasks using the binary segmentation method are applied on sub-samples of the data of differing lengths, and then the results are combined to create an overall answer. We do not restrict the total number of change-points a time series can have, therefore, our proposed method works well when the spacings between change-points are short. Our main change-point detection statistic is the time-varying autoregressive conditional duration model on which we apply a transformation process in order to decorrelate it. To examine the performance of EBS we provide a simulation study for various types of scenarios. A proof of consistency is also provided. Our methodology is implemented in the R package eNchange, available to download from CRAN.


Introduction
Irregularly spaced time series, i.e. data recorded as and when they emerge, have attracted significant attention due to the development and increasing automation in information technology.This type of data, aka high-frequency data, has not only penetrated financial markets, whereby trade transactions, as a result of continuous electronic trading, have generated a vast amount of datasets.Examples of irregularly spaced time series can be found in many industrial and scientific domains such as in natural disasters where floods, earthquakes or volcanic eruptions typically occur at irregular time intervals; in clinical trials where a patient's state of health is typically observed at different points of time; or in data collecting sensors (smart applications of Internet of Things) which are activated only when an event happens to, e.g., preserve battery energy.
Modelling high-frequency data successfully requires that the dynamics in the time series data are captured efficiently and standard time series methods are not appropriate for these data.To address this, a class of models, the autoregressive conditional duration models (henceforth, ACD), was originally introduced by Engle and Russell (1998).In particular, the authors model the conditional mean as an autoregressive process which has a multiplicative relationship with positive valued error terms.These specifications resemble the multiplicative structure of the GARCH model used for modelling an asset's return volatility, but thanks to the positive valued error terms ACDs are more suitable for modelling the dynamics of continuous positive valued random variables, such as durations (the time it takes between two trades occurring in a trading book), trading volume (the size of orders), or market depth (the available liquidity at any point of time).ACD models have become popular with many variations owing partially to the popularity of the GARCH model andits numerous extensions.Without being exhaustive, we refer the reader to Weibull ACD, gamma ACD, and generalized gamma ACD; to models with more flexible dynamic specification such as the logarithmic ACD model of Bauwens and Giot (2000), the Box-Cox ACD model of Dufour and Engle (2000); and to non-parametric or semi-parametric versions found in Cosma and Galli (2006) or Saart et al. (2015).An alternative to ACD's autoregressive specification is to parameterize the intensity function which assumes a self-exciting process.This class of stochastic process was originated by Hawkes (1971) and, hence, these stochastic processes are referred to as Hawkes processes.
They serve as epidemic models: the occurrence of a number of events (e.g.seismic or buy/sell) will increase the probability of further events.Their use in modelling financial high-frequency data has picked up with many applications and variations.
In practice, however, time series entail changes in their dependence structure and the above mentioned models are more appropriate for stationary time series.It will be a crude approximation to adopt them in modeling the non-stationary processes without at least a minimal adaptation in a model's assumptions.There is also high risk in prediction and forecasting in doing so, as noted by Mercurio and Spokoiny (2004) for financial data.Dealing with high-frequency financial data is not less risky.A representative example is the U-shape in a stock's trading activity typically observed during a day.We cannot expect a market's behaviour to be like for like at every point of time.
Information flows almost continuously and, therefore, the intra-day dynamics vary significantly and cannot be ignored.
The simplest type of deviation from stationarity is, arguably, piecewise stationarity.This implies that the parameters of a stochastic process remain constant (hence, the process is stationary) for a certain period of time.In this paper, we focus on this type and we aim to identify the stationary segments by detecting the locations and number of change-points.If, of course, we knew the segments a-priori we would fit a stationary model to each of the segments and proceed with the prediction or forecasting tasks, but this knowledge is typically absent.
Detecting change-points has attracted significant attention.One approach to solve this problem is to formulate it through an optimization task i.e. minimising a multivariate cost function (or criterion) and adding a penalty when the number of change-points is unknown (see (Yao, 1988)).
Depending on the model a user can adopt certain cost functions: the least squares for change in the mean of a series (Yao and Au, 1989) or (Lavielle and Teyssiere, 2007), the Minimum Description Length criterion (MDL) for non-stationary time series (Davis et al., 1995), the Gaussian log-likelihood function for changes in the volatility (Lavielle and Teyssiere, 2007) or the covariance structure of a multivariate time series (Lavielle and Teyssiere, 2006).
Change-point detection methods that adopt a multivariate cost function often come with a high computational cost.Dynamic programming (Bellman and Dreyfus (1966) and Kay (1998)), Segment Neighbourhood (Auger and Lawrence, 1989) or Optimal Partitioning (Jackson et al., 2005), are used in solving change-point problems, but their complexity is at least OpT 2 q where T is the sample size.
An alternative approach to the change-point estimation problem is to minimize a series of univariate cost functions in a 'greedy' manner, i.e. detect a change-point and then progressively move to identify more.A popular representative of this category is the Binary Segmentation method (BS) and the reasons for its popularity are its low computational complexity and easiness of implementation: after identifying a single change-point (through the use of a certain statistic such as the CUSUM) the detection of further change-points continues to the left and to the right of the initial change-point until no further changes are detected.
The BS algorithm has been adopted to solve various types of problems with Inclan and Tiao (1994) to be, perhaps, the first to use it to detect breaks in the variance of a sequence of independent observations.Kim et al. (2000) extend Inclan and Tiao (1994) method to a GARCH(1,1) model and Lee and Park (2001) extend the same method to linear processes.Fryzlewicz and Subba Rao (2011) use the BS method to test for multiple change-points in an ARCH process.Cho and Fryzlewicz (2012) apply the binary segmentation method on the wavelet periodograms in order to identify change-points in the second-order structure of a nonstationary process.Using the wavelet periodogram, Killick et al. (2013) propose a likelihood ratio test under the null and alternative hypotheses, but assume that the number of change-points are bounded from above.The BS method is also used for multivariate (possibly high-dimensional) time series segmentation in Cho and Fryzlewicz (2015) for detecting change-points in the second-order structure and in Cho and Korkas (2018) for detecting change-points in high-dimensional multivariate GARCH processes.For irregularly spaced time series, at least in the context of ACD or Hawkes processes, no literature exists that proposes a BS method (or an alternative change-point detection method) to detect change-points in a piecewise stationary ACD or Hawkes process.We note that Roueff et al. (2016) introduce the time-varying Hawkes process which is locally stationary, not necessarily piecewise stationary and, hence, they do not deal with change-point detection.
The BS algorithm, under specific model specifications, may be inefficient in detecting the change-points.Fryzlewicz (2014) illustrates that with a simple signal+noise model having only three change-points close to each other and in the middle of the series.At the beginning of the BS algorithm the CUSUM statistic does not clearly point to a change-point, hence, it does not move to search for more.This behaviour should not come as a surprise: BS is a greedy algorithm searching for a single change-point at each iteration and failing to do so results in an unnecessary early stop.The Wild Binary Segmentation algorithm (WBS Fryzlewicz ( 2014)), a state-of-theart change-point method, aims to solve the above mentioned limitation using a 'certain random localization mechanism'.We discuss WBS later on, but we note here that this randomized algorithm was extended to univariate time series (Korkas and Fryzlewicz, 2017) and high-dimensional panel data (Wang and Samworth, 2018).
In this work, we capitalize on BS's popularity and propose a new randomised version of it which we term Ensemble Binary Segmentation (EBS).In particular, we draw a number M of random segments from a given univariate time series and apply the BS algorithm in each of these segments.The estimated change-points are then collected from over the M BS applications.Due to this simple mechanism the ways to combine the estimated change-points are numerous, for instance, the final set of change-points comprise of those that appear frequently or appear more frequently relative to other estimated change-points.An extra feature is that the estimated change-points can be ranked based on their relative frequency which is useful when post-processing change-points.
The reader might question why WBS is not preferred for detecting change-points in a piecewise stationary ACD or Hawkes process and a new method is proposed instead.First, EBS exhibits better actual performance both in terms of computational speed and accuracy.The way that WBS works does not allow an efficient post-processing of the estimated change-points resulting in sometimes significant oversegmentation (in the case of a few change-points) or undersegmentation (in the case of frequent change-points).The source of the problem is the spurious detection of change-points which is a consequence of the distributional features of the ACD multiplicative form.Further, EBS, in practice, is fast because it does not search for a single change-point in every iteration, but rather it can locate all the change-points with even a single random draw with a small, albeit non-zero, probability.Second, WBS is a special case of BS while the method we propose here can be extended to include other segmentation algorithms.In other words, when drawing a number of M of random segments we do not have to apply BS in each of these, but any other method can be potentially used.Therefore, studying EBS is a crucial starting point for other ensemble-type of change-point detection techniques.
The action of aggregating many method outputs to randomized versions of the data is not new.
In fact, it has been studied extensively in the statistics and machine learning literature.Random forest is a popular representative that enjoys good predictive performance.Briefly, a random forest works as follows: grow multiple regression (decision) trees to randomized parts of the training dataset, and then average the trees.Our proposed method works in a similar manner and provides a fresh solution to the change-point problem.
The paper is structured as follows: In Section 2 we present the time-varying ACD model (tvACD), the core of our detection algorithm and in Section 3 the time-varying Hawkes process which serves as an alternative to tvACD.The EBS algorithm is presented in Section 4, including its theoretical consistency in detecting the number and locations of change-points.The same section covers the post-processing of the estimated change-points and conducts simulation studies to obtain universal algorithm parameters that almost minimize the need for tuning by the user.Further, in Section 5 we conduct a simulation study to examine the actual performance of our proposed algorithm vis-à-vis the standard BS algorithm.In Section 6 we apply our method to financial highfrequency transaction data.Proofs of the results are in Appendix.Our methodology is implemented in the R package eNchange, available to download from CRAN.
2 Time-varying ACD model Let τ t be the transaction time where 0 " τ 0 ă ... ă τ T ´1 and x t " τ t ´τt´1 be the duration between trades/events.We consider the following time-varying ACD model (1) where ψ t " Erx t |x t , ..., x 1 |θ 1 s is the conditional mean duration of the t-th event with parameter vector θ 1 and Gp.q is a general distribution over p0, `8q with mean equal to 1 and parameter vector θ 2 .In this work we assume that ε t " expp1q, even though extensions to other cases (e.g. the Weibull distribution) are possible but technically challenging.
We denote the vector of parameters involved in the conditional mean ψ t by Θptq " pωptq, α j ptq, β j ptqq P R 1`2j .
Then, Θptq is piecewise constant in t with N change-points η i (η 0 " 0, η N `1 " T ) i.e. at any η i we have that Θpη i q ‰ Θpη i `1q for all t ‰ η i .
The number of change-points N and their locations are assumed to be unknown and we aim to estimate them.The parameter values in Θptq are also unknown, but these can be estimated after the change-points have been detected and stationary segments are identified.
For a non-stationary stochastic process x t , its strong-mixing rate is defined as a sequence of coefficients i.e.
In Fryzlewicz and Subba Rao (2011) the mixing rate of univariate, time-varying ARCH processes was investigated, and in Cho and Korkas (2018) the mixing rate of any pair of time-varying GARCH processes was established.A tvACD process is also strong mixing at a geometric rate, under the following Lipschitz-type condition on the density f εt of ε t .
(A3) The distribution of ε t satisfies the following: for any a ą 0, there exists fixed K ą 0 independent of a such that ż |f εt puq ´fεt pur1 `asq|du ď Ka.
Fryzlewicz and Subba Rao (2011) show that (A3) is satisfied when certain conditions about a density function f hold i.e. the first derivative is bounded, after some finite point m; the derivative f 1 declines monotonically to zero; and ş 8 0 |zf 1 |dz ă 8.It is straightforward to show that f εt satisfies these conditions, hence, (A3) holds when ε t " expp1q.

Alternative model: time-varying Hawkes process
A different approach to modeling the arrival times is through the Hawkes process: a self-exciting process where the arrival of an event (or events) causes the conditional intensity function to increase.

It is defined as follows:
Definifition (Hawkes process) Let tN τ u N ą0 be a counting process on R `with associated history tF τ u τ ě0 .The point process is said to be a Hawkes process if the conditional function λpτ |F τ q takes the following form where λ 0 is the initial intensity at time τ " 0 and V : R Ñ R `is a memory kernel (also termed excitation function).
Kernel Vpτ q modulates the change that event τ t has on the intensity function λpτ q at time τ .
In his original work, Hawkes (1971) assumes an exponential kernel Vpτ q " p ÿ j"1 α j e ´βj τ , α j and β j ą 0 @j. (4) The parameters α j control the increase in the arrival density after each arrival in the system, while parameters β j control the decay rate of this arrival.The unconditional expected value is Epλpτ qq " λ 0 1 ´řp j"1 α j {β j and therefore, in order for stationarity to hold the following condition must be satisfied The integrated intensity function Λ is defined as and from the time-change property the transformed durations τ t ´τt´1 " Λpτ t´1 , τ t q are exponentially distributed with unit rate (or, equivalently, they follow a standard Poisson process).Under the model (3) the integrated intensity can be written as Λpτ t´1 , τ t q " where e ´βj pτt´1τtq for j " 1, . . ., p.
This property is due to the exponential kernel appearing often and, hence, Apt´1q can be recursively computed reducing the calculation of Λpτ t´1 , τ t q (and the likelihood function) from quadratic to linear time.Ogata (1988) calls Λpτ t´1 , τ t q the residual process and it serves as a diagnostic test for the goodness-of-fit of a point process on R `.
We denote the vector of parameters involved in the conditional intensity λpτ q by Θ H pτ q " pλ 0 pτ q, α j pτ q, β j pτ qq P R 1`2j .
Then, Θ H pτ q is piecewise constant in τ with N change-points η i such that 0 ă η 1 ă . . .ă η N ă T (η 0 " 0, η N `1 " T ) i.e. at any η i we have that Θpη i q ‰ Θpη i `1q for all τ ‰ η i .Similar to the tvACD case, we do not know the number of change-points N in Θ H pτ q and their locations which we aim to estimate them.
To prove consistency of EBS in detecting the number of change-points and their locations when the underlying model is a tvHawkes process we would need assumptions around non-negativity of Θ H pτ q and stationarity as in (5) in between any two consecutive change-points.However, to the best of our knowledge, no mixing rates have been established for a time-varying Hawkes process and, therefore, it is not trivial to estimate an upper bound for the cumulative error term.In the next section, we only focus on tvACD while we use tvHawkes as a benchmark (and an alternative) to tvACD.
4 Two-stage change-point detection methodology

Stage A: Transformation of the point processes
In the first stage of our proposed methodology with form the process U t " gpx t , x t´1 , ...x t´p´q q with any fixed p `q.This function g is required to be bounded and Lipschitz continuous.
(A4) The function g : R 1`p`q Ñ R satisfies |g| ď ḡ ă 8 and is Lipschitz continuous, i.e., |gpz 0 , . . ., z p`q q ´gpz 1 0 , . . ., z 1 p`q q| ď C g Empirical residuals are widely adopted for detecting changes in the parameters of a stochastic model and in this work it will form the basic statistic for our detection algorithm in the context of point processes.
For the tvACD model, following Fryzlewicz and Subba Rao (2014) we select g 0 such that where the last term x t is added to ensure the boundness of U t .This transformation decorrelates the original tvACD process and lightens its tails.Therefore, it serves as our main change-point detection statistic by observing that when any parameter Θptq undergoes a change at some point t, so does EpU t q.In this work we favour logpU t ` q where as in ( 7) to reduce the rightward skew observed in the distribution of U t .Our consistency result is based on U t , but the main algorithm uses the log-transformation purely because it gives better results.
It is important to note that any 'diagonal' transformation of x t such as x 2 t or logpx 2 t q should not be used because the squared process x 2 t will be highly autorrelated.This will distort the changepoint detection most likely by producing a false picture of the locations of the change-points.For a discussion on transformations we refer to the original work of Fryzlewicz and Subba Rao (2014).
We prepare the ground for a consistency result for our proposed method.Let tx i t u T t"1 denote a stationary ACD process with parameters Θpη i `1q, and the innovations coinciding with ε t over the associated segment rη i `1, η i`1 s.For each b we also define r U i t " g 0 px i,t´p t , ψ i,t´q t q and denoting the index of the change-point strictly to the left and nearest to t by vptq " maxt0 ď i ď N : η i ă tu with which tx i t u T t"1 and r U i t are defined.
Proposition 2. Suppose that (A2)-(A4) hold, and let y t " U t .Then, we have the following decomposition Proof of Proposition 2 can be found in Appendix.Unlike EpU t q, we have r g t which is exactly constant between any two adjacent change-points without any boundary effects.By its construction, q does not satisfy Epz t q " 0, but thanks to the mixing properties of U t as of Proposition 1, its scaled partial sums can be appropriately bounded.In the next section, we introduce the multiple change-point detection algorithm.

Stage B: The Segmentation algorithm
The Binary Segmentation algorithm We first present the Binary Segmentation algorithm within the framework of the ACD model.In the next section, we explain how this algorithm can be enhanced by ensembling detected change-point from multiple applications of the BS algorithm on smaller (and random) segments.
We consider the CUSUM-type statistic which has been widely adopted for change-point detec- The Ensemble Binary Segmentation algorithm Our recommended enhancement of the BS algorithm, which we argue to be a significant improvement, is based on the fact that the BS method can possibly fit the wrong model when multiple change-points are present as it searches the whole series.The application of the CUSUM statistic can result in spurious change-point detection when e.g. the true change-points occur close to each other.Especially, Olshen et al. (2004) notice that BS method can fail to detect a small change in the middle of a large segment which is illustrated in Fryzlewicz (2014).
To solve this, Fryzlewicz (2014) proposes a randomised version of the binary segmentation method (termed Wild Binary Segmentation -WBS) where the search for change-points proceeds by calculating the CUSUM statistic in smaller segments whose length is random.By doing so, WBS aims to draw favourable intervals with probability tending to one with the sample size containing at most a single change-point.
However, WBS is tailored around BS and the CUSUM statistic and it is not always straightforward to apply to other segmentation methods.In addition, in order to adapt the WBS technique to detecting change-points in, for instance, the second-order structure of a time series as in Korkas and Fryzlewicz (2017), we require a number of new solutions.These include introducing the smallest possible random interval length, and limiting the permitted "unbalancedness" of the CUSUM statistics due to the fact that many of the intervals considered are short, which typically causes spurious behaviour of the corresponding CUSUM statistics.Even with these new solutions in hand, a post-processing step is still needed in order to control the total number of detected change-points.
However, we observe that the post-processing does not contain information about the importance of the estimated change-points.In simple words, a spurious change-point is likely deemed as important as a real change-point, which often resulted in real change-points to be removed after post-processing was applied.We note that this phenomenon does not arise in the signal+iid Gaussian noise setting (Fryzlewicz, 2014), and it is entirely due to the distributional features of the multiplicative setting.
We propose a randomised algorithm, but instead of drawing random intervals, calculate the CUSUM statistic in each of these and proceed in a "to-the-left-and-to-the-right" manner, we run M multiple BS on random segments of the underlying univariate series.We then collect all the obtained change-points and calculate the frequency of the occurrences by simply counting the number of times a certain change-point appears over M draws of the BS algorithm.This results in a better performance compared with BS as it is more likely to draw favourable intervals which can contain a single change-point (as in the WBS methodology) or more than one (as in the BS methodology).By doing so we aim to balance the benefits of both worlds.The action of ensembling (hence, borrowing from the machine learning literature we term our methodology Ensemble BS or EBS) the change-points allows us to rank a change-point based on its importance.
In the last stage, we have the options to either i. inspect the histogram of the estimated change-points; or ii. to keep only the change-points that appear more than rπ z ¨M s times; or iii.
to apply (ii) and then post-process the detected change-points by removing change-points that are ranked lower and/or are 'close' to high ranked change-points.
We chose this setup because change-points in the middle of a long time series is challenging for BS to detect.The CUSUM statistic (in absolute value) will fail to exceed the threshold and the BS algorithm will stop, see Figure 1.On the contrary, when taking random intervals ps m , e m q and then calculating the CUSUM statistic it was more likely for it to exceed the threshold, see the right It is interesting to note that the number of draws M did not alter the change-point performance.
From the bottom two plots in Figure 1 one can see that the shapes of the empirical distributions look identical even though in the first case M " 1000, while in the second M was 10-fold.
We now describe the EBS algorithm in more detail.First, denote by D m " tη  where Ip.q is the indicator function.In the machine learning literature, formula ( 12) is referred as majority voting.We can also obtain the relative frequency of a change-point η i using f ηi " V ηi {|E| to create a relative importance plot (histogram).
Obviously, a change-point that ranks high should be preferred over other change-points.On the other hand, to make a decision about how many change-points are finally selected, one way is to select from the ensembles that η i such that We emphasize that other types of thresholding can be utilized, for instance, a change-point is selected when its rank Rpη i q is above the average rank, i.e.
The relative importance plot also provides a type of 'scree plot', whereby an elbow (kink) indicates what the right number of selected change-points is.This type of scree plot is common in, for example, Principal Component Analysis and the selection of the number of principal components.
It is not hard to see that EBS and WBS are similar in the sense that, for a sufficiently large M , EBS will select a significant number of favourable draws rs m , e m s each containing a single changepoint.Applying the BS method on each of these intervals is equivalent to applying the CUSUM statistic once.In practice, however, we noticed that this is not the case and the added feature of EBS in ranking the estimated change-points resulted in a significant improvement in performance.

Theoretical properties of EBS
In this section we present the consistency theorem for the EBS algorithm for the total number N and locations of the change-points 0 ă η 1 ă ... ă η N ă T ´1 with η 0 " 0 and η N `1 " T .To achieve this, we impose the following conditions in addition to (A1) to (A4) mainly to control the detectability of each η b .
The guaranteed speed of convergence of P pZ T q to 1 is no faster than p1 ´δT T ´1{9q M where M is the number of random draws.
The rate of convergence for the estimated change-points obtained for the BS method by Fryzlewicz and Subba Rao ( 2014) is Op ?T log ᾱ T q where ᾱ ą 0 when δ T is T .In the EBS setting, the rate is square logarithmic when δ T is of order log T , which represents an improvement.
We now elaborate on the minimum number M of random draws required to ensure that the bound on the speed of convergence of P pZ T q to 1 in Theorem 1 is suitably small.Suppose that we wish to ensure that This is equivalent to by noting that logp1 ´yq « ´y around y " 0.
Let us consider the "easiest" case, i.e. δ T " T .This results in a logarithmic number of draws, which leads to particularly low computational complexity, and it also has the same complexity with the WBS case.When δ T " a logpT q, then the required M increases almost linearly, but it has less computational complexity than WBS, which also explains why EBS is generally faster than WBS.
Finally, we discuss π z which appears in ( 13), and it acts as a decision rule when aggregating estimated change-points across M applications of the BS algorithm.Theoretically, π z tends to 0 when T Ñ 8 and for a sufficiently chosen C. In practice, EBS tended to return spurious changepoints due to the distributional features of the multiplicative setting which require us to control the partial sums of z t in (ii) of Proposition 2. Our recommendations for the choice of π z along with the choice of M are discussed in Section 4.2.

Post-processing
In practice, the real number of change-points is not known to us and to reduce the risk of oversegmentation we propose a post-processing method.The need for post-processing the estimated change-points from a detection routine is common within the context of multiplicative models.We refer the reader to Inclan and Tiao (1994), Cho and Fryzlewicz (2012) and Korkas and Fryzlewicz (2017).In these works, the post-processing method compares every detected change-point from the main detection method against the adjacent ones (re)using the CUSUM statistic.Even though there are variations in implementing a post-processing between them, their common drawback is the lack of information about the importance of the detected change-points. is not within distance, we keep this change-point and repeat the process until all change-points are separated by ∆ T .

Choice of parameters for transformation
The right choice of the transformation function g 0 will influence the empirical performance of our methodology and its power in particular.This choice boils down to determine the coefficients C j , j " 0, . . ., p `q.
The BASTA-res algorithm of Fryzlewicz and Subba Rao (2014) performs change-point detection in the univariate ARCH process, as already seen similar to the ACD process, by analysing the transformation of the input time series obtained similarly to U 2 t .They recommend the use of 'dampened' versions of the GARCH parameter estimates which we also adopt here for the ACD process.This leads to the choice of C 0 " p ω, C j " p α j {F, 1 ď j ď p and C p`k " p with within-series dampening parameters F ě 1.
Empirically, the motivation behind the introduction of F is as follows.For x t with time-varying parameters, we often observe that α j and β k are over-estimated in the sense that ř p j"1 p α j `řq is close to one, especially when dealing with real data.There has been evidence in the literature that change-points may cause persistence estimation in volatility models (e.g.Francq et al. (2001)).Mikosch and Stȃricȃ (2004), among others, show that the estimated persistence close to unity in GARCH models is likely spurious and confounded by neglected change-points.We observed the same phenomenon when trying to fit an ACD process to simulated and real data.Using the raw estimates in place of C j 's in (7), therefore, is not the best approach.Cho and Korkas (2018) propose to choose the dampening factor F as By construction, F is bounded as F P r1, 99s and approximately brings p ω and ř p j"1 p α j `řq to the same scale.The selection of the order p, q can be done through the means of an information criterion; however, taking into consideration the possibility that the estimated persistence will be close to unity, we conducted a simulation experiment to establish the right choice of order.
We varied b in βptq above from 0.05 to 0.85, practically meaning a tvACD process with minor to strong persistence in the second segment.For the dampening factor F we chose to vary it from 1 to 10 (instead of using the formula of (Cho and Korkas, 2018)) to understand its impact in changepoint detection performance.For every pair pF, bq we simulated a tvACD process and we performed two types of transformation.In the first case, the unobservable ψ t in the transformation g 0 was replaced with the empirical estimates obtained with the MLE estimates of the ACD parameters.In the second case, we performed the same transformation, but assuming p β 1 " 0 in ( 14).In both cases, we then applied the EBS methodology using the parameters M " 500, π z " 0.05, C 1 " 0.456 (the choice of these default values is discussed in the next section).
The results indicate that the transformation involving the empirical estimates from an ACD(0,1) model vis-a-vis ACD(1,1) had a superior performance.The dampening factor F did not have as a strong impact and generally the performance remained unchanged for increasing F .In a separate experiment not shown here the dynamic selection of F worked better and is, hence, recommended.

Choice of threshold and parameters
In this section we present choices of the parameters involved in the Theorems 1 and 2. In particular, we have that the threshold π T , in both theorems, includes the constant C 1 .To approximate the distribution of the CUSUM statistic in the absence of change-points one approach is to use a parametric resampling procedure that described in Cho and Korkas (2018) and which provides good results.A similar approach has widely been adopted in the change-point literature including Kokoszka and Teyssière (2002) who test for the presence of a change-point in the parameters of univariate GARCH models.
However, this resampling procedure adds computational time and it does not fully take advantage of the dampening transformation aiming to bring a series closer to an iid exponential distribution.For that reason we conduct experiments to establish the universal value of the threshold parameter under the null hypothesis of no change-points.
In particular, we generate stationary ACD processes of size T , varying from 500 to 100000 with a step 50.The exact choice of the ACD model parameters (i.e.pω, α, βq) did not alter the results and, therefore, not reported here.Then we find b that maximises (9).The ratio gives us an insight into the magnitude of parameter C 1 .We repeated this experiment 100 times for different values of T and we selected C 1 as the 100p1 ´αq%-percentile for each instance of T .Our results indicated that C 1 tends to decrease as we increase the sample size and remains unchanged after a certain point (see Figure 2).
To propose a general rule that will apply in most cases we fitted the regression Having estimated the values for ĉ0 , ĉ1 , ĉ2 , ĉ3 we were able to use fitted values for any sample size T .
For samples larger than T " 100000, we used the same C 1 values as for T " 100000.
We turn to the choice of M -the number of ensembles to run -and π z -the relative frequency threshold.In Section 4.2 we discussed that the minimum number of M is in the range of a few thousands when the distance of any two adjacent change-points is of Opplog T q 1{2 q.However, in practice and thanks to the randomised ensemble mechanism of EBS, M can be much smaller.As for π z , the theoretical choice of 0 is adequate, albeit for the distributional features of the multiplicative setting we consider here a non-zero positive value can work better in practice.Similar to the procedure of the choice of C l above, we simulate stationary ACD processes and we apply the EBS algorithm for every triplet pM, π z , T q keeping everything else the same.We vary M from 50 to 5000 with a step 50; π z from 0 to 0.25 with a step 0.01; and T from 5000 to 100000 with a step 1000.Finally, we repeat the experiment, but instead of a stationary process we simulate tvACD processes to establish the right choice of pM, π z , T q in the existence of multiple change-points.In particular, five change-points η 1 , ¨¨¨, η 5 are introduced to the tvACD process parameters in ( 1)-( 2) at 0.25 ¨T, 0.5 ¨T, 0.75 ¨T, 0.85 ¨T, 0.95 ¨T .The parameters change at η 1 , ¨¨¨, η 5 as follows: ωptq " 1.5Ipt ď η 1 q `0.5Ipη 1 ă t ď η 2 q `1.5Ipη 2 ă t ď η 4 q `1Ipt ą η 5 q.
To assess the detection accuracy of our proposed methodology we calculate the ratio of the number of change-points detected within 1% of the real change-points over the number of real change-points.For the model with no change-points, we define detection accuracy as being one when no change-points were detected by EBS; otherwise, we calculate it as 1{p N `1q, where N is the estimated number of change-points returned by EBS.
By inspecting Figure 3, we can see that the choice of π z " 0.05 results in a significant improvement in the detection accuracy compared with lower values.For larger values, we do not see any further improvement, even when the sample size increased.
The most important output of this exercise is that, conditioning on π z , a high number of draws M did not result a considerable improvement in accuracy.On the contrary, a low number in the range of hundreds gave similar results to M " 5000 even for larger samples.
For the simulation study and the real application we set π z " 0.05 and M " 500 which are also the default values in the R package.
5 Simulation study

Models with no change-points
We simulated stationary time series with innovations ε t " expp1q for ACD and sample sizes T " 200, 1000 for different specifications.For the Hawkes process we chose the time horizon h " 500, whereby the sample size will depend on the choice of the parameters.Roughly speaking, for λ 0 " 1, α 1 " 0.1 and β 1 " 0.7 the sample size T « 5000.We report the number of occasions (out of 100) the methods incorrectly rejected the null hypothesis of no change-points for each case.
From Table 5.1, we can see that both BS and EBS performed well meaning that the risk of segmenting a stationary process is limited.It should be mentioned though that when the effective size of a stationary Hawkes process is small and the ratio α 1 {β 1 is approaching 1, both EBS and BS tend to incorrectly reject the null hypothesis of no change-points more often.This was due to the Hawkes estimation optimization routine itself which became erratic and often did not converge.To

Models with change-points
We now examine the detection performance of our method for a set of non-stationary models, both ACD and Hawkes processes.We consider various test models and examine BS and EBS performance over 100 simulations for each of the test model.To compare the performance between the two detection methods we calculate the following error measures: Erp N ´N q 2 s.In addition, the EBS method has improved rates of convergence compared with BS, hence, we examine the total number of change-points identified within t1% ¨T u from the real ones in order to assess the performance on how close the change-points are to their estimates.
Since the purpose of EBS is to also assist in the post-processing stage, its accuracy should be judged in parallel with the total number of change-points identified.We use a test from Korkas and Fryzlewicz (2017) that tries to accomplish this.Assuming that the maximum distance from a real change-point η is denoted by d max , an estimated change-point η is correctly identified if |η ´η| ď d max (here within 1% of the sample size).In the case where two (or more) estimated change-points are within this distance d max then only one change-point, the closest to the real change-point, is classified as correct and the rest are deemed to be false, except if any of these are close to another change-point.An estimator performs well when the hit ratio HR " #correct change-points identified maxpN, N q is close to 1.According to the authors, the term maxpN, N q penalises cases where, for example, the estimator correctly identifies a certain number of change-points all within the distance d max , but N ă N .It also penalises the estimator when overestimates the total number of change-points and all N detected change-points are within distance d max of the true ones.
(M1) tvHawkes processes with a single change-point A change-point η 1 is introduced to the Hawkes process parameters in (6) at 0.95¨h where h " 1000.Parameters λ 0 ptq, α 1 ptq and β 1 ptq Considering this model allows us to assess the detection performance of our proposed methodology when a change-point occurs at the end of a time series which is of interest when data update continuously (on-line).
Both methods performed well by identifying the single change-point within a small distance from the real change-point, and it is worth mentioning that their performance was almost identical.This is an indicative example where selecting EBS over BS is safe even though BS is also expected to be optimal in detecting a single change-point.In this scenario, EBS performed better mainly resulted by the better accuracy in detecting the two change-points.In particular, EBS achieved a hit ratio of 0.93 which was resulted by identifying the two change-points within %1 in %90 cases of the cases.This is a good example where EBS is able to improve accuracy without the risk of oversegmenting a series.
In this setup, EBS performed better than BS again and for the same reason as in M5: similar accuracy between the two methods, but EBS identified the right number of change-points in most cases.
EBS achieved high accuracy here by controlling the total number of the estimated changepoints.We particularly note the reduction in the error measure Erp N ´N q 2 s from 142.85 to 10.81 between BS and EBS, which is also the case for M6.

Application
In this section, we test our proposed methodology using Apple Inc Contract For Differences (CFD) data obtained from Dukascopy for free.CFDs are leveraged products which normally reflect the underlying stock price, but they do not transfer stock ownership to the buyer.They are favoured for their easiness to access (reducing the need to trade in a stock exchange) and the extended trading hours many brokers provide.
The data consists of the bid and ask quotes with their timestamps as published by the broker.
Alternatively, one could choose transaction data which will require further cleansing due to the existence of multiple small size trades.We analyse all the data from June (20 business days) and from 09:30 until 16:00 (US East local time), the stock exchange trading hours.The daily average sample size (number of quotes per business day) is 60917, the minimum 42007, and the maximum 86777.
For each day in month June we ran EBS (assuming a tvACD process) on the raw durations, the time lapsed between two quotes updates.This is to show whether EBS was able to segment a given series without first removing the intraday seasonality, as a piecewise constant model can approximate it well.Our methodology is designed to capture even small deviations from stationarity.
Therefore, it can potentially work well in this setup whereby the underlying dynamics of a duration series change slowly and monotonically at the start of the trading day, flatten during the day and monotonically increase towards the end (resulting in the typical U-shape).
For completeness, we also ran EBS after removing intraday seasonality from the series of durations following Engle and Russell (1998).The seasonality effect is defined as a multiplicative component in the following formula Xt " x t ¨φpτ q where Xt are the raw durations, φpτ q is the intraday seasonality effect and x t are the standardized seasonality durations.To estimate φpτ q we applied a smoothing spline using cross-validation to select the penalty parameter and the degrees of freedom.On the transformed durations x t " Xt { φpτ q, we applied the EBS algorithm.
By inspecting Figures 4 and 5, we see that EBS is able to capture the changes in the trading behaviour during the early and late sessions.Still, when removing seasonality it is very common across the days to observe the existence of change-points between 12.00 and 13.00.We do not argue that these change-points capture seasonality left out from applying a smoothing spline.In fact, we believe that financial high-frequency data are characterized by other various types of nonstationarity which cannot be ignored.This observation has implications both for trading and risk management operations.Given the rise of the electronic trading, intraday risk monitoring will need to adjust in order to incorporate previous change-points and provide more reliable forecasts in anticipation of non-stationarity.On the left, EBS has been applied on the raw durations and on the right the durations has been transformed such that to remove the intraday seasonality.

Conclusion
The work has addressed the problem of detecting the change-points in the structure of an irregularly spaced time series.We proposed a new ensemble-type BS methodology, whereby we combined the change-points detected by applying BS on randomized sub-samples of the data.Doing so, we were able to detect change-points that occur frequently and in short distance between them as well we controlled the total number of change-points that appeared to be spurious.From the simulation study in Section 5 it appears that the EBS mechanism performs well at this task.We believe that, apart from the good performance, EBS is a good starting point of studying other ensemble-type change-point detection methods.where EpV t q ď c 1 ρ t´ηptq´1 1 with ρ 1 P p0, 1q and c 1 ą 0 being some constants independent of t.
Proof.We note that x t being an ACD process does not impact the arguments in the proof of Fryzlewicz and Subba Rao (2014) for an ARCH(q) or Cho and Korkas (2018) for a GARCH(p, q) process, therefore, it is omitted.

Proof of Proposition 2
Part (i) follows trivially from the definition of r U vptq t see Cho and Korkas (2018).
For the proof of (ii), we adopt the arguments similar to those used in the proof of Lemma 4 in Fryzlewicz and Subba Rao (2014).
Recall that vptq denotes the index of the nearest change-point among those satisfying η b ă t, and the definition of z t in the context of time series segmentation, namely z t " U t ´r g t .We want to show that PpF T q Ñ 1, where for some fixed c 1 ą 0.
We investigate the probability of the following event: 1 for λ T " ?log T .For notational convenience let ν " e ´s `1.

Proof of Theorem 1
The proof of consistency is based on the following additive model y t " f t `zt , t " 0, ..., T ´1.
We define the following two CUSUM statistics where n " e ´s `1, the size of the segment defined by ps, eq.

Y b
s,e can be seen as the inner product between sequence ty t u t"s,...,e and a vector ψ b s,e whose elements ψ b s,e,t are constant and positive for t ď b and constant and negative for t ą b such that they sum to zero and sum to one when squared.Similarly for S b s,e .
Let s, e satisfy η p0 ď s ă η p0`1 ă ... ă η p0`q ă e ď η p0`q`1 for 0 ď p 0 ď N ´q.The inequality will hold at all stages of the algorithm until no undetected change-points are remained.We impose at least one of the following conditions s ă η p0`r 1 ´Cδ T ă η p0`r 1 `Cδ T ă e, for some 1 ď r 1 ď q (16) tpη p0`1 ´sq ^ps ´ηp0 qu _ tpη p0`q`1 ´eq ^pe ´ηp0`q qu ď C T where ^and _ denote the minimum and maximum operators, respectively.These inequalities will hold throughout the algorithm until no further change-points are detected.
We In addition, we notice that P `pD M T q c ˘ď P ´pV pM q ηr ě 1q c @ 1 ď r ď N where ζ " 1, ¨¨¨, N K and V pM q ηi as in (12).
For a single (favourable) BS run and on a generic interval ps, eq Ď ps m0 , e m0 q satisfying ( 16) and ( 17 Proof.The case where ps m0 , e m0 q contains a single change-point is shown in Lemma 2 of Korkas and Fryzlewicz (2017) and, therefore, omitted.When two change-points occur then we can show using Lemma 1 from Cho and Fryzlewicz (2012) |S η p 0 `r1 sm 0 ,em 0 | " ˇˇˇˇa η p0`r 1 ´sm0 `1? e m0 ´ηp0`r 1 ?n m0 pf ηi`1 ´fηi q which is bounded from below by C 3 ?δ T when (B3) is satisfied.The above inequality also holds in the case of multiple change-points when a favourable draw n m0 " δ T .
Lemma 3. Assuming that (16) holds, then there exists C 2 ą 0 such that for b satisfying |b ήp0`r 1 | " C 2 γ T for some r 1 , we have |S Proof.This follows the exact same steps as in Fryzlewicz (2014) or Korkas and Fryzlewicz (2017) by taking into consideration Lemma (2) above.
Lemma 4.Under conditions ( 16) and ( 17) there exists 1 ď r 1 ď q such that |b ´ηp0`r 1 | ď T , where b as in (18) and T " C log T for a positive constant C.
Proof.Let f d sm 0 ,em 0 define the best function approximation to f t such that arg max d |xψ d sm 0 ,em 0 , f y| " arg min d ř em 0 t"sm 0 pf t ´f d sm 0 ,em 0 q where f d sm 0 ,em 0 " f `xf, ψ d sm 0 ,em 0 yψ d sm 0 ,em 0 , f is the mean of f and ψ d sm 0 ,em 0 is a set of vectors that are constant and positive until d and then constant and negative from d `1 until e m0 .
If it can be shown that for a certain T ă C 2 γ T , we have em 0 ÿ t"sm 0 pY t ´Ȳ d sm 0 ,em 0 ,t q 2 ą em 0 ÿ t"sm 0 pY t ´f η p 0 `r1 sm 0 ,em 0 ,t q 2 (20) as long as T ď |d ´ηp0`r 1 | then this would prove necessarily that |b ´ηp0`r 1 | ď T .
The rest of the proof follows the same ANOVA-type decomposition from Fryzlewicz (2014).
For the right-hand side of ( 20 For the left-hand side of (20) we apply similar reasoning to the above to bound it i.e.C 5 |d ήp0`r This leads us to the same triplet of inequalities as in the argument in the proof of Theorem 3.2 in Fryzlewicz (2014) i.e. q _ pλ T |d ´ηp0`r 1 | ´1{2 q _ pλ 2 T q. (22) Hence, with the requirement that |d ´ηp0`r 1 | ď C 2 γ T " C 2 λ T ?δ T we obtain and T " maxp1, C 2 qλ T , which concludes the proof.
Lemma 5.Under conditions ( 16) and ( 17) and on the event F T from Proposition 2 where b is given in (18).
Proof.where the third inequality follows from Lemma 3.
Lemma 6.For some positive constants C, C 1 , let s, e satisfy either • D1 ď p ď N such that s ď η p ď e and pη p ´s `1q ^pe ´ηp q ď C T or • D1 ď p ď N such that s ď η p`1 ď e and pη p ´s `1q _ pe ´ηp`1 q ď C 1 T .
Then, on the event F T from Proposition 2 where b is given in (18).
top plot in the same figure.What is important to observe is that only s m has to be as close to the first change-point as possible in order for the CUSUM statistic to surpass the threshold, while end point e m can be much further to the right.This observation holds symmetrically i.e. if e m is close to the last (fourth) change-point then s m can freely take any value in the interval r1, r0.475T sq.In those both favourable scenarios the BS algorithm will proceed to identify the rest change-points consistently.
. .u for m " 1, . . ., M the set of all the change-points detected by the BS algorithm on a random interval rs m , e m s Ď r1, T s.Let E " Ť M m"1 D m , the ensemble collection of all the detected change-points from the M applications of BS, and |E| its cardinality.Evaluate each change-point η i by counting the number of times (frequency) it appears in the ensemble E, i.e.

Figure 1 :
Figure 1: A simulated tvACD process from the model (11) (top left).Various CUSUM statistics calculated on the transformed series U t where the start and end points differ.In blue s m " 1 and e m " 3000; in red s m " 425 and e m P r1515, 2765s; in green s m " 425 and e m P r1375, 2765s (top right).The empirical histogram of the estimated change-points from M " 1000 draws (bottom left).The empirical histogram of the estimated change-points from M " 10000 draws (bottom right).The four vertical bars in red are the real locations of the change-points.

p
ηi for i " 1, . . ., p N Output: p B (B1) The number of change-points N in (1)-(2) is unknown and allowed to increase with T and only the minimum distance between the change-points can restrict the maximum number of N .(B2) There exists a fixed constant f ˚ą 0 such that max 1ďtďT |f t | ď f ˚. (B3) The distance between any two adjacent change-points satisfies min r"1,...,N `1 |η r ´ηr´1 | ě δ T , where δ T ě C log T for a large enough C. (B4) The magnitudes of the change-points satisfy inf 1ďiďN |f ηi`1 ´fηi | ě f ‹ where f ‹ ą 0. Theorem 1 Suppose that Assumptions (A1)-(A5) and (B1)-(B4) hold.With the number of change-points as N and the locations of those change-points as η 1 , ..., η N , let N and η1 , ..., ηN be the number and locations of the change-points (in ascending order) estimated by the Ensemble Binary Segmentation algorithm.There exist constants C 1 and C

Figure 3 :
Figure 3: Heatmaps of detected numbers of change points (measured by the ratio 1{p N `1q) in EBS, depending on M and π z for different sample sizes.A ratio value close to 1 indicates that the exact number of change-points were detected and within 1% ¨T distance from the real ones.

(
M2) tvHawkes processes with two change-points This is a similar model to M1 with the difference that, in addition to the first change-point η 1 that occurs at 0.25 ¨h, a second change-point η 2 occurs at 0.75 ¨h.Further, only the λ 0 ptq parameters changes.3 @ i " 1, ¨¨¨, 3, and β piq 1 " 0.7 @ i " 1, ¨¨¨, 3.The motivation of examining this non-stationary model derives from the seasonality observed in irregularly spaced time series.A typical example is the intraday trading activity of a stock where the majority of trades occurs either in the open or the close sessions.

Figure 4 :
Figure 4: Histogram of detected change-points collected over the sample period grouped by hour.On the left, EBS has been applied on the raw durations and on the right the durations has been transformed such that to remove the intraday seasonality.

Figure 5 :
Figure 5: Sample days of Apple Inc CFD durations before adjusted for seasonality (left) and after seasonality is removed (right).In red, the estimated change-points using EBS.
first define N K ď N clusters of the N change-points C ζ where 1 ď ζ ď N such that in each cluster the distance between the first (from the left) and the last change-point is of c ¨δT where c ą 0. We note that a cluster can contain a single change-point.In the case of multiple change-points where the distance between any two adjacent change-points is larger than δ T there will be N K clusters with |C 1 | " ¨¨¨" |C N | " 1.We now define intervals (not necessarily symmetric) I L ζ and I R ζ around clusters such that for every triplet tC ζ´1 , C ζ , C ζ`1 u ˙ for ζ " 1, ..., N K where δ ζ min " min tmintC ζ u ´maxtC ζ´1 u, mintC ζ`1 u ´maxtC ζ uu.We recall that when applying the EBS algorithm M intervals ps m , e m q, m " 1, ..., M are drawn from a discrete uniform distribution over r1, T s.These M intervals also define M runs of the BS algorithm each of which can contain none, a single or multiple change-points.We define the event D M T as D M T " t@ζ " 1, ..., N K D m " 1, ..., M ps m , e m q P I L ζ ˆIR ζ u.
y t ¸for s ď b ă e. (9) A large value of |Y b s,e | typically indicates the presence of a change-point in the level of y t in the vicinity of t " b.In particular, if Y s,e ą π T where Y s,e " max and π T is a threshold (the choice of which is discussed in Section 4.2), then the location of the Algorithm 1: BinSeg (Binary Segmentation algorithm) Input: ty t u, π T , s, e, p B Step 1: compute Y b s,e for s ď b ă e; Step 2: set Y b s,e Ð max sďbăe |Y b s,e | and p η Ð arg max sďbăe Y b s,eStep 3: if Y s,e ą π T then p B Ð p B Y tp ηu p B Ð BinSeg(ty t u, π T , s, p η, p B) p B Ð BinSeg(ty t u, π T , p η `1, e, p B) end Output: p B

Table 1 :
Stationary processes results.For each of the Hawkes processes the average sample size T a vg is given in brackets.Figures show the number of occasions the methods incorrectly detected change-points.

Table 2 :
Nonstationary processes change-point detection results.The error measures are selfexplanatory while the hit ratio HR is described in the main text.