Modeling and forecasting time series is of fundamental importance in various applications. There may be occasional changes of behavior in a time series. Some examples are the changes in the stock market due to the financial crisis, or the variations of an EEG signal caused by the mode change in the brain. In the econometrics literature, this kind of time series is referred to as regime-switching model [1, 2]. In regime switching models, the time series is assumed to have states, and if belongs to state (
), the probability density function (pdf) ofconditioning on its past is in the form of . The autoregressive (AR) model, one of the commonly used techniques to model stationary time series , is usually used to model each state. The autoregression of state is given by where
are independent and identically distributed (i.i.d.) noises with zero mean and variance. Here , ,
is a real-valued vector of lengththat characterizes state . A more detailed survey on this model can be found in . We refer to this model as a multi-state AR model and to as the AR filter or AR coefficients of state . The above model with was first analyzed by Lindgren  and Baum et al. . The model with general is widely studied in the speech recognition literature 
. The multi-state AR model is a general statistical model that can be used to fit data in many real world applications. It was shown that the model is capable of representing non-linear and non-stationary time series with multimodal conditional distributions and with heteroscedasticity
. There are two basic underlying assumptions in this model: 1. Autoregression assumption, which is reasonable if the observations are obtained sequentially in time; 2. Multi-state assumption, which is reasonable if the stochastic process exhibits different behaviors in different time epochs. For example, stock prices may have dramatic while not permanent changes in the case of business cycles or financial crises, and those dynamics can be described by stochastic transitions among different states.
Despite the wide applications of the multi-state AR model, there are few results on how to estimate the number of statesin a time series. Obviously, different values of produce a nested family of models and models with larger ’s fit the observed data better. The drawback of using complex models with a large
is the over-fitting problem which decreases the predictive power of the model. Hence, a proper model selection procedure that identifies the appropriate number of states is vital. It is tempting to test the null hypothesis that there arestates against the alternative of . Unfortunately, the likelihood ratio test of this hypothesis fails to satisfy the usual regularity conditions since some parameters of the model are unidentified under the null hypothesis. An alternative is to apply Akaike information criterion (AIC)  or Bayesian information criterion (BIC)  to introduce a penalty on the complexity of the model in the model selection procedure. However, in general AIC and BIC are shown to be inaccurate in estimating the number of states .
In this paper, we propose a model selection criterion inspired by the work of Tibshirani et al.  who studied the clustering of i.i.d. points under Euclidean distance. The idea is to identify by comparing the goodness of fit for the observed data with its expected value under a null reference distribution. To that end, we first draw a reference curve which plots the “goodness of fit” versus based on the most non-informative distributed data, and describes how much adding new AR states improves the goodness of fit. We then draw a similar curve based on the observed data. In this work we choose the “goodness of fit” measure to be the mean squared prediction error (MSPE). Finally, the point at which the gap between the two curves is maximized is chosen as the estimated .
Besides the simplicity and effectiveness, another benefit of the proposed model selection criterion is that it is adaptive to the underlying characteristics of AR processes. The criterion for the processes of little dependency, i.e., the roots of whose characteristic polynomial are small, is different from the criterion for those of large dependency. In this sense, it takes into account the characteristics behind the observed data in an unsupervised manner, even though no domain knowledge or prior information is given.
formulates a specific class of the multi-state AR model, where the transitions between the states are assumed to be a first order Markov process. We emphasize that this parametric model is considered primarily for simplicity and the proposed Gap statistics can be applied to general multi-state AR processes. A new initialization approach is also proposed that can effectively reduce the impact of a bad initialization on the performance of the expectation-maximization (EM) algorithm. SectionIV presents some numerical results to evaluate the performance of the proposed approach. Experiments show that the accuracy of the proposed approach in estimating the number of AR states surpasses those of AIC and BIC.
Ii Gap Statistics
This section describes our proposed criterion for selecting the number of states in a multi AR process, inspired by . We draw a reference curve, which is the expected value of MSPE under a null reference distribution versus , and use its difference with the MSPE of the observed data to identify the number of states, . We show that computing each point of the reference curve turns out to be a clustering problem in the space of AR coefficients of a fixed size, where the distance measure for clustering is derived from the increase in MSPE when a wrong model is specified. We derive the distance measure in closed form, introduce an approach to generate stable AR filters that are uniformly distributed, and apply the -medoids algorithm to approximate the optimal solution for the clustering problem. We first outline our proposed model selection criterion in Subsection II-A, and then elaborate on the distance measure in Subsections II-B and the generation of random AR filters in Subsections II-C.
Ii-a The Model Selection Criterion
We use superscript to represent the data at time step , and
to denote the normal distribution with meanand variance . Symbols in bold face represent vectors or matrices. We start from a simple scenario where the data is generated using a single stable AR filter : , where , , , and are i.i.d. . Suppose we are at time step and we want to predict the value at time . If is used for prediction, the MSPE is . But if another AR filter is used for prediction instead of , i.e., , the MSPE becomes . The difference of the two MSPE is defined by
It is easy to observe that is always nonnegative, which means that using the mismatch filter for prediction increases MSPE. We refer to as the mismatch distance between two filters and , though it is not a metric. When the data generated from has zero mean, i.e., , we let also represents of length (with constant term omitted) with a slight abuse of notation, and we use in the same manner.
As has been mentioned in Section I, our model selection criterion is based on a reference curve that describes how much adding a new state increases the goodness of fit in the most non-informative or the “worst” case. To that end, we consider an -state zero mean AR process where at each time step , nature chooses random mismatch filters (with zero constants) for prediction. In such a worst scenario, the filters that minimize the average mismatch distances to the random filters are naturally believed to be the true data generating filters, and that minimal value, which is the average MSPE, is plotted as the reference curve. This leads to the following clustering problem in the space of stable AR filters , where
Clustering of Stable Filters: For a fixed , let , , , be a set of uniformly generated stable filters of a given length . We cluster into disjoint clusters , and define the within cluster sum of distances to be
where is defined in (1) and will be further simplified in (4), (5) and (6). By computing for , we obtain the reference curve. The optimization problem (2) can be solved by the -medoids algorithm .
The model selection criterion is outlined in Table 1. We note that the bound for the roots is determined by the estimated filters, and thus the reference is data-dependent. Intuitively, if the process has less dependency, or in other words a point has less influence on its future points, the roots of the characteristic polynomials of each AR process are closer to zero and the MSPE curve will have smaller values. Thus, the filters from which the reference curve is calculated should also be drawn from a smaller bounded space.
Ii-B Distance Measure for Autoregressive Processes
In this subsection, we provide the explicit formula for the distance in Equation (1). Assume that the data is generated by a stable filter of length . Let be the characteristic polynomial of , and let denote the roots of , i.e., , where lie inside the unit circle (). Similarly define for . The value in (1) can be computed using the power spectral density and Cauchy’s integral theorem as:
for , where denotes the complex conjugate of . For the degenerate cases when or , reduces to or .
For now we assume that at each state has zero mean by default, unless explicitly pointed out. We use in Identity (4) instead of in Identity (3) to compute the reference curve. The derived reference curve can be applied to the general case. The reason is that it is more difficult to detect two AR states with the same mean than those that have different means. Therefore, the reference curves for the zero mean case (the “worst” case) can be used in general.
The distance between two AR filters can be explicitly expressed in terms of the coefficients. This is computationally desirable if the filters are random samples generated in the coefficient domain, as will be discussed in Subsection II-C.
Notations: Consider two polynomials of nonnegative powers and respectively of degrees and . Let respectively denote the reciprocal polynomial of , and the multiplication of and , i.e., , . Let be the resultant of and . Define and , where are the roots of .
The values of and can be computed as polynomials of the coefficients of and .
The proof follows from the fact that the resultant of and is given by the determinant of their associated Sylvester matrix , and that for any , can be computed as polynomials in the coefficients of via Newton’s identities. We further provide the following result.
Another simple way to compute the distance measure is given by the following lemma.
Let be the true filter of an autoregression with zero mean. The variance , the correlations , and the covariance matrix of the autoregression are respectively defined to be Define , for and , if and otherwise . Then and can be computed by
where is determined by . The value of in terms of and can be computed by
Ii-C Generating Uniformly Distributed Filters with Bounded Roots
As mentioned before, Gap statistics requires a reference curve that is calculated by clustering the filters randomly chosen from a reference distribution. In some scenarios we need to generate sample filters from , where is calculated from the observed data. Inspired by the work of Beadle and Djurić , we provide the following result on how to generate a random point in with uniform distribution
Generation of an independent uniform sample of can be achieved by the following procedure:
1. Draw uniformly on the interval ;
2. For , suppose that we have obtained that is uniformly distributed in . Draw independently from a pdf proportional to the following function on the interval
We prove by induction. The pdf of is proportional to one. For , suppose that the pdf of is proportional one inside and zero elsewhere. Suppose that and are determined by (8). The Levinson-Durbin recursion in (8) automatically enforces the stability constraint that falls inside . The pdf of can be computed as
where is the Jacobian from to () taking to be given. Therefore, if is proportional to the value given by (7), the joint pdf of is proportional to one in and zero elsewhere. ∎
A sample of that is uniformly distributed in can be generated by the recursion , where and are independently generated.
A popular way to describe the switching behavior between different states is to assume that the transition between the states follows a first-order Markov process. In this section, we adopt this assumption to formulate a parametric multi-state AR model for illustration purpose, even though the model selection criterion proposed in Section II is applicable to other multi-state AR models.
Iii-a Notations and Formulations
Let denote the set of data points that are generated from state . Suppose that are fixed and known. Let and be a sequence of missing (unobserved) indicators, where is a matrix, is a vector, and
Clearly, . We note that is a binary vector of length containing a unique “”; with a slight abuse of notation is the location of that “”. We assume that
is a Markov chain with transition probability matrix, where , and is drawn from , where denotes the family of multinomial distributions. In other words, the assumed data generating process (given a fixed ) is:
Let be the set of unknown parameters to be estimated, where is of length (including the constant term). Though computing the maximum-likelihood estimation (MLE) of the above probabilistic model (10) is not tractable, it can be approximated by a local maximum via the EM algorithm . The EM algorithm produces a sequence of estimates by the recursive application of E-step and M-step to the complete log-likelihood until a predefined convergence criterion is achieved. The complete log-likelihood can be written as
For brevity, we provide the EM formulas below without derivation. In the E-step, we obtain a function of unknown parameters by taking the expectation of (11) with respect to the missing data and given the most updated parameters,
can be computed recursively. We note that the parameters involved in the right-hand side of (13) take values from the last update. In the M-step, we use the coordinate ascent algorithm to obtain the following local maximum. The “old” superscriptions are omitted for brevity.
Iii-B Initialization of EM
The convergence speed of the EM algorithm strongly depends on the initialization and an improper initialization can cause it to converge to a local maximum which is far away from the global optimum. A routine technique is to use multiple random initializations and choose the output with the largest likelihood , but this can be significantly time consuming. Here, we use a new initialization technique to get a fast and reliable convergence for the EM algorithm. This technique is based on the fact that for time series obtained in most practical areas, the self-transition probability of the states is usually close to one, i.e., . By adopting this assumption, we propose the initialization method in Algorithm 4, which is shown empirically to produce more reliable and efficient EM results. We note that the “split” style rule that appears in line 5 of Algorithm 4 is used elsewhere (e.g. s).