Covariance-based Dissimilarity Measures Applied to Clustering Wide-sense Stationary Ergodic Processes
Qidi Peng Nan Rao Ran Zhao
stat.ML 62-07, 60G10, 62M10







































































Covariance-based Dissimilarity Measures Applied to Clustering Wide-sense Stationary Ergodic Processes
    Qidi Peng[Institute of Mathematical Sciences, Claremont Graduate University, Claremont, CA 91711. Email: qidi.peng@cgu.edu.] Nan Rao[School of Mathematical Sciences, Shanghai Jiao Tong University, Shanghai, China. Email: nan.rao@sjtu.edu.cn.] Ran Zhao[Institute of Mathematical Sceinces and Drucker School of Management, Claremont Graduate University, Claremont, CA 91711. Email: ran.zhao@cgu.edu.]
    
==================================================================================================================================================================================================================================================================================================================================================================================================================




We introduce a new unsupervised learning problem: clustering wide-sense stationary ergodic stochastic processes. A covariance-based dissimilarity measure together with asymptotically consistent algorithms is designed for clustering offline and online datasets, respectively. We also suggest a formal criterion on the efficiency of dissimilarity measures, and discuss of some approach to improve the efficiency of our clustering algorithms, when they are applied to cluster particular type of processes, such as self-similar processes with wide-sense  stationary ergodic increments. Clustering synthetic data and real-world data are provided as examples of applications.



Keywords:  cluster analysis · wide-sense stationary ergodic processes · covariance-based dissimilarity measure · self-similar processes

MCS (2010):  62-07 ·  60G10 · 62M10

















§ INTRODUCTION


Cluster analysis, as a core category of unsupervised learning techniques, allows to discover hidden patterns in data where one does not know the true answer upfront.  Its goal is to assign a heterogeneous set of objects into non-overlapping clusters, where in each cluster any two objects are more related to each other than to objects in other clusters. Given its exploratory nature, clustering has nowadays a number of applications in various fields of both industry and scientific research, such as biological and medical research <cit.>, information technology <cit.>,  signal and image processing <cit.>, geology <cit.> and finance <cit.>. There exists a rich literature of cluster analysis on random vectors, where the objects, waiting to be clustered, are sampled from high-dimensional joint distributions. There is no shortage of such clustering algorithms <cit.>. However, stochastic processes are quite a different setting from random vectors, since their observations (sample paths) are sampled from processes distributions. While the cluster analysis on random vectors has developed aggressively, clustering on stochastic processes receives much less attention. Today cluster analysis on stochastic processes deserves increasingly intense study, thanks to their vital importance to many applied areas, where the collected information are indexed by real time and are especially long. Examples of these time-indexed information include biological data, financial data, marketing data, surface weather data, geological data and video/audio data, etc. 

Recall that in the setting of random vectors, a process of clustering often consists of two steps:


  Step 1 One suggests a suitable dissimilarity measure to describe the distance between 2 objects, under which two objects are close to each other becomes meaningful.

  Step 2 One designs an enough accurate and computationally efficient clustering function based on the above dissimilarity measure.

Clustering stochastic processes is performed in a similar way but new challenges may arise in both Step 1 and Step 2. Intuitively, one can always apply existing  random vectors clustering approaches to cluster arbitrary stochastic processes, such as non-hierarchical approaches (K-means clustering methods) and hierarchical approaches (agglomerative method, divisive method) <cit.>, based on naive dissimilarity measures (e.g., Euclidean distance, Manhattan distance or Minkowski distance). However, one faces at least 2 potential risks when applying the above approaches to clustering stochastic processes:


  Risk 1 These approaches might suffer from their huge complexity costs, due to the great length of their sample paths. As a result classical clustering algorithms are often computationally forbidding <cit.>.

  Risk 2 These approaches might suffer from over-fitting issues. For example, clustering stationary or periodic processes based on Euclidean distance between the paths, without considering their path properties will result in over fitting, bad clusters situation.

In summary, classical dissimilarity measures or clustering strategies would fail in clustering stochastic processes.

Fortunately, the complexity cost and the over-fitting errors of clustering processes could be largely reduced, if one is aware of the fact that a stochastic process often possesses fine paths features (e.g., stationarity, Markov property, self-similarity, sparsity, seasonality, etc.), which is unlike an arbitrary random vector. An appropriate dissimilarity measure then should be chosen to be able to capture these paths features. Clustering processes is then performed to group any two sample paths into one group, if they are relatively close to each other under that particular dissimilarity measure.  Below are some examples provided in the literature.

<cit.> proposed a dissimilarity measure between two special sample paths of processes. In their setting it is supposed that, for each path only sparse and irregularly spaced measurements with additional measurement errors are available. Such features occur commonly in longitudinal studies and online trading data. Based on this particular dissimilarity measure, classification and cluster analysis could be made. <cit.> developed a new algorithm to perform clustering of multivariate and functional data, based on a covariance-based dissimilarity measure. Their attention is focused on the specific case of a set of observations from two populations, whose probability distributions have equal mean but differ in terms of covariances. <cit.> designed consistent algorithms for clustering strict-sense stationary ergodic processes (see the forthcoming Eq. (<ref>) for the definition of strict-sense ergodicity), where the dissimilarity measure is proposed as distance of process distributions. It is worth noting that the consistency of their algorithms is guaranteed thanks to the assumption of strict-sense ergodicity.

In this framework, we aim to design asymptotically consistent algorithms to cluster a general class of stochastic processes, i.e., wide-sense stationary ergodic processes (see Definition <ref> below). Asymptotically consistent algorithms can be obtained for this setting, since the covariance stationarity and ergodicity allow the process to present some featured asymptotic behavior with respect to their length, rather than to the total number of paths. 


A stochastic process X={X_t}_t∈ T (the time indexes set T can be either ℝ_+=[0,+∞) or ℕ={1,2…}) is called wide-sense stationary if its mean and covariance structure are finite and time-invariant: 𝔼(X_t)=μ for any t∈ T, and for any subset (X_i_1,…,X_i_r), its covariance matrix remains invariant subject to any time shift h>0:

    ℂov(X_i_1,…,X_i_r)= ℂov(X_i_1+h,…,X_i_r+h).

Denote by γ the auto-covariance function of X. Then X is further called weakly ergodic (or wide-sense ergodic) if it is ergodic for the mean and the second-order moment:


  * If X is a continuous-time process (e.g., T=ℝ_+), then  it satisfies for any s∈ℝ_+, 

    1/h∫_s^s+hX_u uμ,

and

    1/h∫_s^s+h(X_u+τ-μ)(X_u-μ) uγ(τ), ,

where  denotes the almost sure convergence (convergence with probability 1).

  * If X is a discrete-time process (e.g., T=ℕ), then  it satisfies for any s∈ℕ∪{0}, 

    X_s+X_s+1+…+X_s+h/h+1μ,

and

    ∑_u=s^s+h(X_u+τ-μ)(X_u-μ)/h+1γ(τ), .



Wide-sense stationarity and ergodicity are believed to be a very general assumption, at least in the following senses:


  * The assumption that each process is generated by some mean and covariance structure is sufficient for capturing all features of a wide-sense stationary ergodic process. In other words, our algorithms intend to cluster means and auto-covariance functions, not process distributions. 

  * Wide-sense stationary ergodic process partially extends the strict-sense one. A finite-variance strict-sense stationary ergodic process (see Eq. (<ref>) for its definition) is also wide-sense stationary ergodic. However strict-sense stationary ergodic stable processes are not wide-sense stationary, because their variances explode <cit.>.

  * A Gaussian process can be fully identified only by its mean and covariance structure. Then a wide-sense stationary  ergodic Gaussian process is also strict-sense stationary ergodic.

  * In the clustering problem, the dependency among the sample paths can be arbitrary. 

There is a long list of processes which are wide-sense stationary ergodic, but not necessarily stationary in the strict sense. The examples of wide-sense stationary processes below are not exhausted.


  Example 1 Non-independent White Noise.

Let U be a random variable uniformly distributed over (0,2π) and define

    Z(t):=√(2)cos(tU),  t∈ℕ.

The process Z={Z(t)}_t∈ℕ is then a white noise because it verifies 

    𝔼(Z(t))=0, 𝕍ar(Z(t))=1  ℂov(Z(s),Z(t))=0,  s≠ t.

We claim that Z is wide-sense stationary ergodic, which can be obtained by using the Kolmogorov's strong law of large numbers, see e.g. Theorem 2.3.10 in <cit.>. However Z is not strict-sense stationary since

    (Z(1),Z(2))≠ (Z(2),Z(3)) .

Indeed, it is easy to see that

    0<𝔼(Z(1)^2Z(2))≠𝔼(Z(2)^2Z(3))=0.


  Example 2 Auto-regressive Models.

It is well-known that an auto-regressive model {Y(t)}_t∼ AR(1) in the form:

    Y(t)=aY(t-1)+Z(t), |a|<1,a≠ 0,  t∈ℕ

is wide-sense stationary ergodic. However it is not necessarily strict-sense stationary ergodic, when the joint distributions of the white noise {Z(t)}_t are not invariant with time-shifting (e.g., take {Z(t)}_t to be the white noise in Example 1).

  Example 3 Increment Process of Fractional Brownian Motion.

Let {B^H(t)}_t be a fractional Brownian motion with Hurst index H∈(0,1) (see <cit.>). For each h>0, its increment process {Z^h(t):=B^H(t+h)-B^H(t)}_t is finite-variance strict-sense stationary ergodic <cit.>. As a result it is also wide-sense stationary ergodic. More detail will be discussed in Section <ref>.


  Example 4 Increment Process of More General Gaussian Processes.

<cit.> introduced a general class of zero-mean Gaussian processes X={X(t)}_t∈ℝ having stationary increments. Its variogram ν(t):=2^-1𝔼(X(t)^2) satisfies:


  (1) There is a non-negative integer d such that ν is 2d-times continuously differentiable over [-2,2], but not 2(d+1)-times continuously differentiable over [-2,2].

  (2) There are 2 real numbers c≠ 0 and s_0∈(0,2), such that for all t∈[-2,2],

    ν(t)=ν^(2d)(0)+c|t|^s_0+r(t),

where the remainder r(t) satisfies:


  * r(t)=o(|t|^s_0), as t→0.

  * There are two real numbers c' > 0, ω > s_0 and an integer q > ω+1/2 such that r is q-times continuously differentiable on [-2, 2]\{0} and for all t ∈ [-2, 2]\{0}, we have

    |r^(q)(t)|≤ c'|t|^ω-q.



It is shown that the process X extends fractional Brownian motion and it also has wide-sense (and strict-sense) stationary ergodic increments when d+s_0/2∈(0,1) (see Proposition 3.1 in <cit.>).

The problem of clustering processes via their means and covariance structures leads us to formulating our clustering targets in the following way.


Let 

    G={G_1,…,G_κ}

be a partitioning of ℕ={1,2,…} into κ disjoint sets G_k, k=1,…,κ, such that the means and covariance structures of 𝐱_i, i∈ℕ are identical, if and only if i∈ G_k for some k=1,…,κ. Such G is called ground-truth of covariance structures. We also denote by G|_N the restriction of G to the first N sequences:

    G|_N={G_k∩{1,…,N}: k=1,…,κ}.


Our clustering algorithms will aim to output the ground-truth partitioning G, as the sample length grows. Before stating these algorithms, we introduce the inspiring framework done by <cit.>.



 §.§ Preliminary Results: Clustering Strict-sense Stationary Ergodic Processes


<cit.> considered the problem of clustering strict-sense stationary ergodic processes. The main fruit in <cit.> is obtaining the so-called asymptotically consistent algorithms to cluster processes of that type. We briefly state their work below. Depending on how the information is collected, the stochastic processes clustering problems consist of dealing with two models: offline setting and online setting.


  Offline setting: The observations are assumed to be a finite number N of paths: 

    𝐱_1 = (X_1^(1),…, X_n_1^(1)),…,𝐱_N = (X_1^(N),…, X_n_N^(N)).

Each path is generated by one of the κ different unknown process distributions. In this case, an asymptotically consistent clustering function should satisfy the following.


A clustering function f is consistent for a set of sequences S if f(S,κ)=G. Moreover, denoting n=min{n_1,…,n_N}, f is called strongly asymptotically consistent in the offline sense if with probability 1 from some n on it is consistent on the set S, i.e.,

    ℙ(lim_n→∞f(S,κ)=G)=1.

It is called weakly asymptotically consistent if lim_n→∞ℙ(f(S,κ)=G)=1.


  Online setting: In this setting the observations, having growing length and number of scenarios with respect to time t, are denoted by 

    𝐱_1 = (X_1^(1),…, X_n_1^(1)),…,𝐱_N(t) = (X_1^(N(t)),…, X_n_N(t)^(N(t))),

where the index function N(t) is non-decreasing with respect to t. 

Then an asymptotically consistent online clustering function is defined below:


A clustering function is strongly (RESP. weakly) asymptotically consistent in the online sense, if for every N∈ℕ the clustering f(S(t),κ)|_N is strongly (RESP. weakly) asymptotically consistent in the offline sense, where f(S(t),κ)|_N is the clustering f(S(t),κ) restricted to the first N sequences:

    f(S(t),κ)|_N={f(S(t),κ)∩{1,…,N}: k=1,…,κ}.



There is a detailed discussion on the comparison of offline and online settings in <cit.>, stating that these two settings have significant differences, since using the offline algorithm in the online setting by simply applying it to the entire data observed at every time step, does not result in an asymptotically  consistent algorithm. Therefore separately and independently studying these two settings becomes necessary and meaningful. 

As the main results in <cit.>, asymptotically consistent clustering algorithms for both offline and online settings are designed. They are then successfully applied to clustering synthetic and real data sets.

Note that in the framework of <cit.>, a key step is introduction to the so-called distributional distance <cit.>: the distributional distance between a pair of process distributions ρ_1, ρ_2 is defined to be

    d(ρ_1,ρ_2)=∑_m,l=1^∞w_m w_l∑_B∈ B^m,l|ρ_1(B)-ρ_2(B)|,

where:


  * The sets B^m,l, m,l≥1 are obtained via the partitioning of ℝ^m into cubes of dimension m and volume 2^-ml, starting at the origin.

  * The sequence of weights {w_j}_j≥1 is positive and decreasing to zero. Moreover it should be chosen such that the series in (<ref>) is convergent. The weights are often suggested to give precedence to earlier clusterings, protecting the clustering decisions from the presence of the newly observed sample paths, whose corresponding distance estimates may not yet be accurate. For instance, it is set to be w_j=1/j(j+1) in <cit.>.

Further, the distance between two sample paths 𝐱_1, 𝐱_2 of stochastic processes is given by

    d(𝐱_1,𝐱_2)=∑_m=1^m_n∑_l=1^l_nw_m w_l∑_B∈ B^m,l|ν(𝐱_1,B)-ν(𝐱_2,B)|,

where:


  * m_n,l_n (≤ n) can be arbitrary sequences of positive integers  increasing to infinity, as n→∞. 

  * For a process path  𝐱=(X_1,…,X_n), and an event B, ν(𝐱,B) denotes the average times that the event B occurs over n-m_n+1 time intervals. More precisely,

    ν(𝐱,B):=1/n-m_n+1∑_i=1^n-m_n+1 1{(X_i,…,X_i+m_n-1)∈ B}.


The process distribution X from which 𝐱 is sampled is called strictly ergodic if

    ℙ(lim_n→∞ν(𝐱,B)=ℙ(X∈ B))=1, .

The assumption that the processes are ergodic leads to that d is a strongly consistent estimator of d: 

    ℙ(lim_n→∞d(𝐱_1,𝐱_2)=d(ρ_1,ρ_2))=1,

where ρ_1,ρ_2 are the process distributions corresponding to 𝐱_1,𝐱_2, respectively.

Based on the distances d and their estimates d, the asymptotically consistent algorithms for clustering stationary ergodic processes in each of the offline and online settings are provided (see Algorithms 1, 2 and Theorems 11, 12 in <cit.>). <cit.> also show that their methods can be implemented efficiently: they are at most quadratic in each of their arguments, and are linear (up to log terms) in some formulations.



 §.§ Statistical Setting: Clustering Wide-sense Stationary Ergodic Processes

Inspired by the framework of <cit.>, we consider the problem of clustering wide-sense stationary ergodic processes.  We first introduce the following covariance-based dissimilarity measure, which is one of the main contributions of this paper. 
(Covariance-based dissimilarity measure)

The covariance-based dissimilarity measure d^* between a pair of processes X^(1), X^(2) (in fact X^(1), X^(2) denote two covariance structures, each may contain different process distributions) is defined as follows:

    d^*(X^(1),X^(2)):= ∑_m,l = 1^∞ w_m w_l
          ×ℳ(
    (𝔼(X_l… l+m-1^(1)),
    ℂov(X_l… l+m-1^(1))), 
    (
    𝔼(X_l… l+m-1^(2)), ℂov(X_l… l+m-1^(2)))),

where:


  * For j=1,2, {X_l^(j)}_l∈ℕ denotes some path sampled from the process X^(j). We assume that all possible observations of the process X^(j) is a subset of {X_l^(j)}_l∈ℕ. For l'≥ l≥ 1, we define the shortcut notation  X_l… l'^(j):=(X_l^(j),…,X_l'^(j)).

  * The function ℳ is defined by: for any p_1,p_2,p_3∈ℕ, any 2 vectors v_1,v_2∈ℝ^p_1 and any 2 matrices A_1,A_2∈ℝ^p_2× p_3,

    ℳ((v_1,A_1),(v_2,A_2)):= |v_1-v_2|+ρ^*(A_1,A_2).


  * The distance ρ^* between 2 equal-sized matrices M_1,M_2 is defined to be

    ρ^*(M_1,M_2):=M_1-M_2_F,

with ·_F being the Frobenius norm:

for an arbitrary matrix M={M_ij}_i=1,…,m; j=1,…,n,

    M_F:=√(∑_i=1^m∑_j=1^n|M_ij|^2).

Introduction to the matrices distance ρ^* is inspired by <cit.>. The matrices distance given in <cit.> is used to measure the distance between 2 correlation matrices. However, our distance ρ^* is a modification of the one in the latter paper. Indeed, unlike <cit.>, ρ^* is a well-defined metric distance, as it satisfies the triangle inequalities. 

  * The sequence of positive weights {w_j} is chosen such that d^*(X^(1),X^(2)) is finite. Observe that the distances |·| and ρ^* in (<ref>)  do not depend on l, as a result we necessarily have

    ∑_l=1^∞ w_l<+∞.

In practice a typical choice of weights we suggest is w_j=1/j(j+1), j=1,2,…. This is because, for most of the well-known covariance stationary ergodic processes (causal ARMA(p,q), increments of fractional Brownian motions, etc.), their auto-covariance functions are absolutely summable: denote by γ_X the auto-covariance function of {X_t}_t,

    ∑_h=-∞^+∞|γ_X(h)|<+∞.

<cit.> pointed out that (<ref>) is a sufficient condition for {X_t} being mean-ergodic. However (<ref>) does not necessarily imply that {X_t} is covariance-ergodic. It becomes a sufficient and necessary condition if {X_t} is Gaussian. 
Therefore subject to (<ref>), taking w_j=1/j(j+1), we obtain for any integer N>0,

    ∑_m,l = 1^N w_m w_l |𝔼(X_l… l+m-1^(1))-𝔼(X_l… l+m-1^(2))|
          = ∑_m,l = 1^N w_m w_l √(m)|μ_1-μ_2|=|μ_1-μ_2|∑_m,l = 1^N1/√(m)(m+1)l(l+1)
          ≤ |μ_1-μ_2|∑_m,l = 1^+∞1/√(m)(m+1)l(l+1)<+∞,

with μ_j=𝔼(X_1^(j)), for j=1,2; and

    ∑_m,l = 1^N w_m w_l ρ^*(ℂov(X_l… l+m-1^(1)),ℂov(X_l… l+m-1^(2)))
          ≤∑_m,l = 1^N w_m w_l √(2∑_k_1=1^m∑_k_2=1^m(γ_X(|k_1-k_2|))^2)
          = ∑_m,l = 1^N w_m w_l √(2∑_q=-(m-1)^m-1(m-|q|)(γ_X(|q|))^2)
          ≤∑_m,l = 1^N w_m w_l √(2m∑_q=-(m-1)^m-1(γ_X(|q|))^2)
          ≤ c∑_m,l = 1^N√(2m)/m(m+1)l(l+1)≤ c∑_m,l = 1^+∞√(2m)/m(m+1)l(l+1)<+∞,

where the constant c=∑_q=-∞^∞|ℂov(X_1,X_1+|q|)|<+∞. Therefore combining (<ref>) and (<ref>) leads to

    d^*(X^(1),X^(2))<+∞.

Hence d^*(X^(1),X^(2)) in (<ref>) is well-defined.


In (<ref>) and (<ref>) we see that the behavior of the dissimilarity measure d^* is jointly explained by the Euclidean distance of means and the matrices distance of covariance matrices. If the means of the processes X^(1) and X^(2) are priorly known to be equal, the distance d^* can be simplified to:

    d^*(X^(1),X^(2))= ∑_m,l = 1^∞ w_m w_lρ^*(
    ℂov(X_l… l+m-1^(1)), ℂov(X_l… l+m-1^(2))).

Note that this dissimilarity measure can be applied on self-similar processes, since they are all zero-mean (see Section <ref>). 

Next we provide consistent estimator of d^*(X^(1),X^(2)).
For 1≤ l≤ n and m≤ n-l+1, define μ^*(X_l… n, m) to be the empirical mean of a process X's sample path (X_l,…,X_n):  

    μ^*(X_l… n,m) :=1/n-m-l+2∑_i = l ^n-m+1 (X_i … X_i+m-1)^T,

and define ν^*(X_l… n, m) to be the empirical covariance matrix of (X_l,…,X_n):  

    ν^*(X_l… n,m)    :=   1/n-m-l+2∑_i = l ^n-m+1 (X_i … X_i+m-1)^T(X_i … X_i+m-1)
          -μ^*(X_l… n,m)μ^*(X_l… n,m)^T,

where M^T denotes the transpose of the matrix M.

Recall that the notion of wide-sense ergodicity is given in Definition <ref>. The ergodicity theorem concerns what information can be derived from an average over time about the ensemble average at each point of time. For the wide-sense stationary ergodic process X, being either  continuous-time or discrete-time, the following statement holds: every empirical mean μ^*(X_l… n,m) is a strongly consistent estimator of the path mean 𝔼(X_l… l+m-1); and every empirical covariance matrix ν^*(X_l… n,m) is a strongly consistent estimator of the covariance matrix ℂov(X_l… l+m-1) under the Frobenius norm, i.e., for all m≥1, we have

    ℙ(lim_n→∞|μ^*(X_l… n,m)-𝔼(X_l… l+m-1)|=0)=1

and

    ℙ(lim_n→∞ν^*(X_l… n,m)-ℂov(X_l… l+m-1)_F=0)=1.

Next we introduce the empirical covariance-based dissimilarity measure d^*, serving as a consistent estimator of the covariance-based dissimilarity measure d^*.

Given two processes' sample paths 𝐱_j=(X_1^(j),…,X_n_j^(j)), j=1,2. Let n=min{n_1,n_2}, we define the empirical covariance-based dissimilarity measure between 𝐱_1 and 𝐱_2  by

    d^*(𝐱_1,𝐱_2):=∑_m= 1^m_n∑_l= 1^n-m+1 w_m w_l
          ×ℳ( (μ^*(X^(1)_l… n,m), ν^*(X^(1)_l… n,m)),(μ^*(X^(2)_l… n,m), ν^*(X^(2)_l… n,m))).

The empirical covariance-based dissimilarity measure between a sample path 𝐱_i and a process X^(j) (i,j∈{1,2}) is defined by

    d^*(𝐱_i,X^(j)):=∑_m= 1^m_n∑_l= 1^n-m+1 w_m w_l
          ×ℳ( (μ^*(X^(i)_l… n,m), ν^*(X^(i)_l… n,m)),(𝔼(X_l… l+m-1^(j)), ℂov(X_l… l+m-1^(j)))).


Unlike the dissimilarity measure d^* which describes some distance between stochastic processes, the empirical covriance-based dissimilarity measure is some distance between two sample paths (finite-length vectors). We will show in the forthcoming Lemma <ref> that d^* is a consistent estimator of d^*.

Two observed sample paths possibly have distinct lengths n_1, n_2, therefore in  (<ref>) we consider computing the distances between their subsequences of length n=min{n_1,n_2}. In practice we usually take m_n=⌊log n⌋, the floor number of log n. 

It is easy to verify that both d^* and d^* satisfy the triangle inequalities, thanks to the fact that both the Euclidean distance and  ρ^* satisfy the triangle inequalities. More precisely, the following holds.


Thanks to (<ref>) and the definitions of d^* (see (<ref>)) and d^* (see (<ref>)), we see that the triangle inequality holds for the covariance-based dissimilarity measure d^*, as well as for its empirical estimate d^*. Therefore for arbitrary processes X^(i), i = 1,2,3 and arbitrary finite-length sample paths 𝐱_i, i=1,2,3, we have

    d^*(X^(1),X^(2))  ≤ d^*(X^(1),X^(3)) + d^*(X^(2),X^(3)),
    
          d^*(𝐱_1,𝐱_2)  ≤d^*(𝐱_1,𝐱_3) + d^*(𝐱_2,𝐱_3),
    
          d^*(𝐱_1,X^(1)) ≤d^*(𝐱_1,X^(2)) + d^*(X^(1),X^(2)).


Remark <ref> together with the fact that the processes are weakly ergodic, leads to Lemma <ref> below, which is the key to demonstrate that our clustering algorithms in the forthcoming section are asymptotically consistent.
 
Given two paths

    𝐱_1=(X_1^(1),…,X_n_1^(1))     𝐱_2=(X_1^(2),…,X_n_2^(2)),

sampled from the wide-sense stationary ergodic processes X^(1) and X^(2) respectively, we have

    ℙ(lim_n_1,n_2 →∞d^*(𝐱_1,𝐱_2) = d^*(X^(1),X^(2)))=1

and

    ℙ(lim_n_i →∞d^*(𝐱_i,X^(j)) = d^*(X^(1),X^(2)))=1,  i,j∈{1,2}, i≠ j.



We take n=min{n_1,n_2}. To show (<ref>) holds it suffices to prove that for arbitrary ε>0, there is an integer N>0 such that for any n≥ N, with probability 1,

    |d^*(𝐱_1,𝐱_2) - d^*(X^(1),X^(2))|<ε.

Define the sets of indexes

    S_1(n)={(m,l)∈ℕ^2: m≤ m_n, l≤ n-m+1}  S_2(n)=ℕ^2\ S_1(n).

To be more convenient we also denote by

    V(X_l… l+m-1^(j)):=(𝔼(X_l… l+m-1^(j)),ℂov(X_l… l+m-1^(j)))

and

    V(X_l… n^(j),m):=(μ^*(X_l… n^(j),m),ν^*(X_l… n^(j),m)),

for (m,l)∈ℕ^2 and j=1,2. By using the definitions of d^* (see (<ref>)), of d^* (see (<ref>)) and the triangle inequality 
    |∑_i∈ Ia_i|≤∑_i∈ I|a_i|, ,
 we obtain

    | d^*(𝐱_1, 𝐱_2) - d^*(X^(1), X^(2))|
          = | ∑_(m,l)∈ S_1(n) w_m w_l (ℳ ( V(X_l… n^(1),m),V(X_l… n^(2),m) ).
          - ∑_(m,l)∈ S_1(n)∪ S_2(n) w_m w_lℳ (V(X_l… l+m-1^(1)), V(X_l… l+m-1^(2)) ) | 
          ≤| ∑_(m,l)∈ S_1(n) w_m w_l (ℳ ( V(X_l… n^(1),m),V(X_l… n^(2),m) ).
          .- ℳ (V(X_l… l+m-1^(1)), V(X_l… l+m-1^(2)) )) | 
          +∑_(m,l)∈ S_2(n)w_mw_lℳ (V(X_l… l+m-1^(1)), V(X_l… l+m-1^(2)) ) 
          ≤∑_(m,l)∈ S_1(n) w_m w_l |ℳ ( V(X_l… n^(1),m),V(X_l… n^(2),m) )
          - ℳ (V(X_l… l+m-1^(1)), V(X_l… l+m-1^(2)) )| 
          +∑_(m,l)∈ S_2(n)w_mw_lℳ (V(X_l… l+m-1^(1)), V(X_l… l+m-1^(2))).

Next note that the metric ℳ satisfies the following triangle inequality:

    |ℳ ( V(X_l… n^(1),m),V(X_l… n^(2),m) )- ℳ (V(X_l… l+m-1^(1)), V(X_l… l+m-1^(2)) )|
          ≤ℳ ( V(X_l… n^(1),m),V(X_l… l+m-1^(1)) )+ ℳ (V(X_l… n^(2),m), V(X_l… l+m-1^(2))).

It follows from (<ref>) and (<ref>) that

    | d^*(𝐱_1, 𝐱_2) - d^*(X^(1), X^(2))|
          ≤∑_(m,l)∈ S_1(n) w_m w_l (ℳ ( V(X_l… n^(1),m),V(X_l… l+m-1^(1)) )
          + ℳ (V(X_l… n^(2),m), V(X_l… l+m-1^(2)) )) 
          +∑_(m,l)∈ S_2(n)w_mw_lℳ (V(X_l… l+m-1^(1)), V(X_l… l+m-1^(2)) ).

Next we show that the right-hand side of (<ref>) converges to 0 as n→∞. First observe that the weights {w_m}_m≥1 have been chosen such that

    ∑_m,l = 1^∞ w_m w_l ℳ(V(X_l… l+m-1^(1)),V(X_l… l+m-1^(2)))<+∞.

Then for arbitrary fixed ε>0, we can find an index J such that for n≥ J,

    ∑_(m,l)∈ S_2(n) w_mw_lℳ(V(X_l… l+m-1^(1)),V(X_l… l+m-1^(2)))≤ε/3.

Next, the weak ergodicity of the processes X^(1) and X^(2) implies that: for each (m,l)∈ℕ^2, V(X_l… n^(j),m) (j=1,2) is a strongly consistent estimator of  V(X_l… l+m-1^(j)), under the metric ℳ, i.e., with probability 1,

    lim_n→∞ℳ(V(X_l… n^(j),m),  V(X_l… l+m-1^(j)))=0.

Thanks to (<ref>), for any (m,l)∈ S_1(J), there exists some N_m,l (which depends on m,l) such that for all n≥ N_m,l, we have, with probability 1, 

    ℳ ( V(X^(j)_l… n,m), V(X_l… l+m-1^(j)) )≤ε/3w_m w_l#S_1(J),  j = 1, 2,

where # A denotes the number of elements included in the set A. 
Denote by N_J=max_(m,l)∈ S_1(J)N_m,l. Then observe that, for n≥max{N_J,J},

    ∑_(m,l)∈ S_2(n)w_mw_lℳ (V(X_l… l+m-1^(1)), V(X_l… l+m-1^(2)) )
          ≤∑_(m,l)∈ S_2(J)w_mw_lℳ (V(X_l… l+m-1^(1)), V(X_l… l+m-1^(2)) ).

It results from (<ref>), (<ref>), (<ref>) and (<ref>) that, for n≥max{N_J,J},

    | d^*(𝐱_1, 𝐱_2) - d^*(X^(1), X^(2))|
          ≤∑_(m,l)∈ S_1(n) w_m w_l  ℳ(V(X^(1)_l… n,m),V(X_l… l+m-1^(1)) )
          +∑_(m,l)∈ S_1(n) w_m w_l ℳ ( V(X^(2)_l… n,m),V(X_l… l+m-1^(2)))
          +∑_(m,l)∈ S_2(J)w_mw_lℳ (V(X_l… l+m-1^(1)), V(X_l… l+m-1^(2)) )
          ≤ε/3+ε/3+ε/3=ε,

which proves (<ref>). The statement (<ref>) can be proved analogously.





§ ASYMPTOTICALLY CONSISTENT CLUSTERING ALGORITHMS
 


 §.§ Offline and Online Algorithms

In this section we introduce the asymptotically consistent algorithms for clustering offline and online datasets respectively. We explain how the two algorithms work, and prove that both algorithms are asymptotically consistent. It is worth noting that the asymptotic consistency of our algorithms relies on the assumption that  the number of clusters κ is priorly known. The case for κ being unknown has been studied in <cit.> in the problem of clustering strictly stationary ergodic processes. However in the setting of wide-sense stationary ergodic processes, this problem remains open.

Algorithm <ref> below presents the pseudo-code for clustering offline datasets. It is a centroid-based clustering approach. One of its main features is that the farthest 2-point initialization applies. The algorithm selects the first two cluster centers by picking the two farthest observations among all observations (Lines 1 - 3), under the empirical dissimilarity measure d^*. Then each next cluster center is chosen to be the observation farthest to all the previously assigned cluster centers (Lines 4 - 6). Finally the algorithm assigns each remaining observation to its nearest cluster (Lines 7-11).



We point out that Algorithm <ref> is different from Algorithm 1 in <cit.> at two points:


  * As mentioned previously, our algorithm relies on the covariance-based dissimilarity d^*, in lieu of the process distributional distances.

  * Our algorithm suggests 2-point initialization, while Algorithm 1 in <cit.> randomly picks 1-point as the first cluster center. The latter initialization was proposed for use with k-means clustering by <cit.>. Algorithm 1 in <cit.> requires κ N distance calculations, while our algorithm requires N(N-1)/2 distances calculations. It is very important to point out that, to reduce the computational complexity cost of our algorithm, it is fine to replace our 2-point initialization with the one in <cit.>. However there are two reasons based on which we recommend using our approach of initialization:


  Reason 1 In the forthcoming Section <ref>, our empirical comparison to <cit.> shows that the 2-point initialization turns out to be more accurate in clustering than the 1-point initialization.

  Reason 2 Concerning the complexity cost, we have the following loss and earn: on one hand, the 2-point initialization requires more steps of calculations than the 1-point initialization; on the other hand, in our covariance-based dissimilarity measure d^* defined in (<ref>), the matrices distance ρ^* requires m_n^2 computations of Euclidean distances, while the distance  ∑_B∈ B^m,l|ν(𝐱_1,B)-ν(𝐱_2,B)| 
given in (<ref>) requires at least n_1+n_2-2m_n+2 computations of Euclidean distances (see Eq. (33) in <cit.>). Note that we take m_n=⌊log n⌋ (⌊·⌋ denotes the floor integer number) though this framework. Therefore the computational complexity of the covariance-based dissimilarity d^* makes the overall complexity of Algorithm <ref> quite competitive to the algorithm in <cit.>, especially when the paths lengths n_i, i=1,…,n are relatively large, or when the database of all distance values are at hand.


Next we present the clustering algorithm for online setting. As mentioned in <cit.>, one regards recently-observed paths as unreliable observations, for which sufficient information has not yet been collected, and for which the estimators of the covariance-based dissimilarity measures are not accurate enough. Consequently, farthest-point initialization would not work in this case; and clustering on all available data results in not only mis-clustering unreliable paths, but also in clustering incorrectly those for which sufficient data are already available. The strategy is presented in Algorithm <ref> below: clustering based on a weighted combination of several clusterings, each obtained by running the offline algorithm (Algorithm <ref>) on different portions of data.

More precisely, Algorithm <ref> works as follows. Suppose the number of clusters κ is known. At time t, a sample S(t) is observed (Lines 1 - 2), the algorithm iterates over j= κ,…,N(t) where at each iteration Algorithm <ref> is utilized to cluster the first j paths in S(t) into κ clusters (Lines 6 - 7). For each cluster its center is selected as the observation having the smallest index among that cluster, and their indexes are ordered increasingly (Line 8). The minimum inter-cluster distance γ_j (see <cit.>) is calculated as the minimum distance d^* between the κ cluster centers obtained at iteration j (Line 9).  Finally, every observation in S(t) is assigned to the nearest cluster, based on the weighted combination of the distances between this observation and the candidate cluster centers obtained at each iteration on j (Lines 14 - 17).


In Algorithm <ref>, β(j) denotes a function indexed by j, which is the value chosen for the weight w_j. Remark that for online setting, our algorithm requires the same number of distance calculations as in Algorithm 2 in <cit.>. They are both bounded by 𝒪(N(t)^2). Using 2-point initialization, our Algorithm <ref> then takes advantage in the overall computational complexity cost. Finally we note that both Algorithm <ref> and Algorithm <ref> require κ≥2. When κ is known, this restriction is not a practical issue. 



 §.§ Consistency and Computational Complexity of the Algorithms

In this section we prove the asymptotic consistency of Algorithms <ref> and <ref>.  They are stated in the 2 theorems below.


Algorithm <ref> is strongly asymptotically consistent (in the offline sense), provided that the true number κ of clusters is known, and each sequence 𝐱_i, i = 1,…,N is sampled from some wide-sense stationary ergodic process.



Similar to the idea used in the proof of Theorem 11 in <cit.>, to prove the consistency statement we will need Lemma <ref> to show that if the sample paths in S are long enough, the sample paths that are generated by the same process covariance structure are closer to each other than to the rest. Therefore, the sample paths chosen as cluster centers are each generated by a different covariance structure, and since the algorithm assigns the rest to the closest clusters, the statement follows. More formally, let n_min denote the shortest path length in S:

    n_min:= min{n_i: i=1,…,N}.

Denote by δ_min the minimum non-zero covariance-based dissimilarity measure between any 2 covariance structures:

    δ_min:= min{d^*(X^(k),X^(k')): k,k'∈{1,…,κ}, k≠ k'}.

Fix ε∈ (0, δ_min/4). Since there are a finite number N of observations, by Lemma <ref> there is n_0 such that for n_min≥ n_0 we have

    max_l ∈{1,…,κ}
    i ∈ G_l ∩{ 1,…,N }d^*(𝐱_i, X^(l))≤ε,

where G_l, l = 1,…,κ denote the covariance structure ground-truth partitions given by Definition <ref>. 

On one hand, by using (<ref>), the triangle inequality (see Remark <ref>) and the fact that 

    max_i∈ I(a_i+b_i)≤max_i∈ Ia_i+max_i∈ Ib_i

for any indexes set I and any real numbers a_i's and b_i's, we obtain

    max_l ∈{1,…,κ}
    i,j ∈ G_l ∩{ 1,…,N }d^*(𝐱_i, 𝐱_j)
          ≤max_ l ∈{1,…,κ}
    i,j ∈ G_l ∩{ 1,…,N }d^*(𝐱_i, X^(l))+max_l ∈{1,…,κ}
    i,j ∈ G_l ∩{ 1,…,N }d^*(𝐱_j, X^(l))
          = max_ l ∈{1,…,κ}
    i ∈ G_l ∩{ 1,…,N }d^*(𝐱_i, X^(l))+max_l ∈{1,…,κ}
    j ∈ G_l ∩{ 1,…,N }d^*(𝐱_j, X^(l))
          ≤ 2 ε< δ_min/2.

On the other hand, by using the triangle inequality (see Remark <ref>), (<ref>) and (<ref>), we have for n_min≥ n_0,

    min_ k,k'∈{1,…,κ},k≠ k'
     i ∈ G_k ∩{ 1,…,N }
    j ∈ G_k'∩{ 1,…,N }d^*(𝐱_i, 𝐱_j) 
          ≥min_k,k'∈{1,…,κ},k≠ k'
     i ∈ G_k ∩{ 1,…,N }
    j ∈ G_k'∩{ 1,…,N }{ d^*(X^(k), X^(k')) - d^*(𝐱_i, X^(k))- d^*(𝐱_j, X^(k'))}
          ≥δ_min-2ε> δ_min/2.

In words, (<ref>) together with (<ref>) indicates that the sample paths in S that are generated by the same covariance structure are closer to each other than to the rest of sample paths. Then by (<ref>) and (<ref>), for n_min≥ n_0, we necessarily have each sample path should be close enough to its cluster center, i.e.,

    max_i = 1,…,Nmin_k = 1,…,κ-1d^* (𝐱_i, 𝐱_c_k)>δ_min/2,

where the κ cluster centers' indexes c_1,…,c_κ are given by Algorithm <ref> as

    (c_1,c_2) := _i,j = 1,…,N, i<jd^* (𝐱_i, 𝐱_j),

and

    c_k :=_i = 1,…,Nmin_j = 1,…,k-1d^* (𝐱_i, 𝐱_c_j), k = 3,…,κ.

Hence, the indexes c_1,…, c_κ will be chosen to index the sample paths generated by different process covariance structures. Then by (<ref>) and (<ref>), each remaining sample path will be assigned to the cluster center corresponding to the sample path generated by the same process covariance structure. Finally Theorem <ref> results from (<ref>), (<ref>) and (<ref>).




Algorithm <ref> is strongly asymptotically consistent (in the online sense), provided the true number of clusters κ is known, and each sequence 𝐱_i, i ∈ℕ is sampled from some wide-sense stationary ergodic process.



The idea of the proof is similar to that of Theorem 12 in <cit.>. The main  differences between the 2 proofs are made by the fact that our covariance-based dissimilarity measure d^* is not bounded by some constant. Although it is not mentioned in the pseudo-code Algorithm <ref>, the notations γ_j's and η are dependent of t, therefore we denote γ_j^t:=γ_j and η^t:=η through this proof.  
In the first step, by using the triangle inequality we can show that for any t>0, any N∈ℕ,

    sup_j∈{1,…,N}
     k∈{1,…,κ}d^*(𝐱_j^t ,X^(k))≤sup_j∈{1,…,N}
     k∈{1,…,κ}(d^*(X^(k), X^(k_j')) +d^*(𝐱_j^t, X^(k_j')) )
          ≤sup_j∈{1,…,N}
     k∈{1,…,κ}d^*(X^(k), X^(k_j')) +sup_j∈{1,…,N}
     k∈{1,…,κ}d^*(𝐱_j^t, X^(k_j'))
          = sup_j∈{1,…,N}
     k∈{1,…,κ}d^*(X^(k), X^(k_j')) +sup_j∈{1,…,N}d^*(𝐱_j^t, X^(k_j')),

where for each j, k_j' is chosen such that 𝐱_j^t is sampled from the process covariance structure X^(k_j'). On one hand,
let

    δ_max:= max{d^*(X^(k),X^(k')): k,k'∈{1,…,κ}, k≠ k'},

then the first term on the right-hand side of (<ref>) can be bounded by the constant δ_max, which neither depends on t nor on N:

    sup_j∈{1,…,N}
     k∈{1,…,κ}d^*(X^(k), X^(k_j')) ≤δ_max.

On the other hand, since 𝐱_j^t is sampled from X^(k_j'), by using the weak ergodicity (see Lemma <ref>), for j=1,…,N, with probability 1,

    lim_t→∞d^*(𝐱_j^t, X^(k_j'))=0.

This together with the fact that a convergent sequence is also bounded, leads to, for each j∈{1,…, N}, there is b_j (not depending on t) such that

    d^*(𝐱_j^t, X^(k_j'))≤ b_j, .

Therefore the second term on the right-hand side of (<ref>) can be bounded as:

    sup_j∈{1,…,N}d^*(𝐱_j^t, X^(k_j'))≤max{b_1,…,b_N}.

Let 

    B(N):=δ_max+max{b_1,…,b_N}.

It is important to point out that B(N) depends only on N but not on t. It follows from (<ref>),  (<ref>), (<ref>) and (<ref>) that

    sup_j∈{1,…,N}
     k∈{1,…,κ}d^*(𝐱_j^t ,X^(k))≤ B(N).

Let δ_min be the one given in (<ref>). Fix ε∈ (0, δ_min/4). By using (<ref>), we can choose some J>0 so that 

    ∑_j = J+1^∞ w_j≤ε.

Recall that in online setting, the ith sample path's length n_i(t) grows with time, for each i. Therefore, by the wide-sense ergodicity (see Lemma <ref>), for every j ∈{1,…,J} there exists some T_1(j)>0 such that for all t ≥ T_1(j) we have

    max_k ∈{1,…,κ}
     i ∈ G_k ∩{ 1,…,j }d^*(𝐱_i^t, X^(k)) ≤ε.

For k = 1,…,κ, define s_k(N(t)) to be the index of the first path in S(t) sampled from the covariance structure X^(k), i.e.,

    s_k(N(t)) := min{ i ∈ G_k ∩{1,…,N(t)}}.

Note that s_k(N(t)) depends only on N(t). 
Then denote

    m(N(t)) := max_k ∈{1,…,κ} s_k(N(t)).

By Theorem <ref> for every j ∈{m(N(t)),…,J} there exists some T_2(j) such that (S(t)|_j, κ) is asymptotically consistent for all t ≥ T_2(j), where S(t)|_j = {𝐱_1^t, …, 𝐱_j^t } denotes the subset of S(t) consisting of the first j sample paths.  Let

    T:= max_i=1,2
     j ∈{1,…,J} T_i(j).

Recall that, by the definition of m(N(t)) in (<ref>), S(t)|_m(N(t)) contains sample paths from all κ distinct covariance structures. Therefore, similar to obtaining (<ref>), for all t ≥ T, we use the triangle inequality, (<ref>) and (<ref>) to obtain

    min_ k,k'∈{1,…,κ}
     k≠ k'd^*(𝐱_c_k^m(N(t))^t , 𝐱_c_k'^m(N(t))^t )
          ≥min_ k,k'∈{1,…,κ}
     k≠ k'(d^*(X^(k), X^(k')) -(d^*(𝐱_c_k^m(N(t))^t, X^(k)) + d^*(𝐱_c_k'^m(N(t))^t, X^(k'))  ))
          ≥δ_min - 2ε≥δ_min/2.

From Algorithm <ref> (see Lines 9, 11) we see 

    η^t : = ∑_j=1^N(t)w_jγ_j^t,    γ_j^t:=min_k,k'∈{1,…,κ}
     k≠ k'd^*(𝐱_c_k^j^t,𝐱_c_k'^j^t).

Hence, by (<ref>), for all t≥ T,

    η^t ≥w_m(N(t))δ_min/2.

For j∈{ J+1,…,N(t)}, by the triangle inequality and (<ref>), we have for all t≥ T,

    γ_j^t=min_ k,k'∈{1,…,κ}
     k≠ k'd^*(𝐱_c_k^j^t , 𝐱_c_k'^j^t )
          ≤min_ k,k'∈{1,…,κ}
     k≠ k'(d^*(X^(k), X^(k')) +(d^*(𝐱_c_k^j^t, X^(k)) + d^*(𝐱_c_k'^j^t, X^(k'))  ))
          ≤δ_max + 2B(N(t)).

Denote by

    M(N(t)):=δ_max + 2B(N(t)),

then (<ref>) can be interpreted as: for all t≥ T, 

    sup_j∈{J+1,…,N(t)}γ_j^t≤ M(N(t)).

By (<ref>), (<ref>) and (<ref>), for every k ∈{1,…,κ} we obtain

    1/η^t∑_j=1^N(t) w_j γ_j^td^*(𝐱_c_k^j^t ,X^(k))
          =1/η^t∑_j=1^J w_j γ_j^td^*(𝐱_c_k^j^t ,X^(k))+1/η^t∑_j=J+1^N(t) w_j γ_j^td^*(𝐱_c_k^j^t ,X^(k))
          ≤1/η^t∑_j=1^J w_j γ_j^t d^*(𝐱_c_k^j^t ,X^(k)) + 2B(N(t))M(N(t))/w_m(N(t))δ_min∑_j=J+1^N(t)w_j
          = 1/η^t∑_j=1^m(N(t))-1 w_j γ_j^t d^*(𝐱_c_k^j^t ,X^(k)) +1/η^t∑_j=m(N(t))^J w_j γ_j^t d^*(𝐱_c_k^j^t ,X^(k))
          + 2B(N(t))M(N(t))ε/w_m(N(t))δ_min.

Next we provide upper bounds of the first 2 items in the right-hand side of (<ref>). On one hand, by the definition of m(N(t), the sample paths in S(t)|_j for j = 1,…,m(N(t)) - 1 are generated by at most κ -1 out of the κ process covariance structures. Therefore for each j ∈{1,…,m(N(t))-1} there exists at least one pair of distinct cluster centers that are generated by the same process covariance structure. Consequently, by (<ref>) and the definition of η^t, for all t ≥ T and k∈{1,…,κ},

    1/η^t∑_j=1^m(N(t))-1 w_j γ_j^t d^*(𝐱_c_k^j^t ,X^(k)) ≤ε/η^t∑_j=1^m(N(t))-1 w_j γ_j^t ≤ε.

On the other hand, since the clusters are ordered in the order of appearance of the distinct covariance structures, we have 𝐱_c_l^j^t = 𝐱_s_l(N(t))^t for all j = m,…,J and l = 1,…,κ, where the index s_l(N(t)) is defined in (<ref>). Therefore, by (<ref>) and the definition of η^t, for all t ≥ T and every l=1,…,κ we have

    1/η^t∑_j=m(N(t))^Jw_j γ_j^t d^*(𝐱_c_l^j^t,X^(l)) =  d^*(𝐱_s_l(N(t))^t ,X^(l)) 1/η^t∑_j=m(N(t))^J w_j γ_j^t ≤ε.

Combining (<ref>), (<ref>), (<ref>) and (<ref>) we obtain, for t≥ T,

    1/η^t∑_j=1^N(t) w_j γ_j^t d^*(𝐱_c_k^j^t ,X^(k))≤ε(2 + 2B(N(t))M(N(t))/w_m(N(t))δ_min)

for all l = 1,…,κ.

Now we explain how to use (<ref>) to prove the asymptotic consistency of Algorithm <ref>. Consider an index i ∈ G_k' for some k' ∈{1,…,κ}. Then on one hand, using (<ref>) and (<ref>), we get for k∈{1,…,κ}, k≠ k',

    1/η^t∑_j=1^N(t) w_j γ_j^t d^*(𝐱_i^t,𝐱_c_k^j^t)
          ≥1/η^t∑_j=1^N(t) w_j γ_j^t d^*(𝐱_i^t ,X^(k)) - 1/η^t∑_j=1^N(t) w_j γ_j^t d^*(𝐱_c_k^j^t ,X^(k)) 
          ≥1/η^t∑_j=1^N(t) w_j γ_j^t (d^*(X^(k),X^(k'))- d^*(𝐱_i^t, X^(k')))
          - 1/η^t∑_j=1^N(t) w_j γ_j^t d^*(𝐱_c_k^j^t ,X^(k)) 
          ≥δ_min - 2ε(2 + 2B(N(t))M(N(t))/w_m(N(t))δ_min).

On the other hand, for any N∈ℕ, by using the wide-sense ergodicity, there is T(N) such that for all t≥ T(N),

    max_k ∈{1,…,κ}
     i ∈ G_k ∩{ 1,…,N }d^*(𝐱_i^t, X^(k)) ≤ε.

Since ε can be arbitrarily chosen, it follows from (<ref>) and (<ref>) that

    _k ∈{1,…,κ}1/η^t∑_j=1^N(t) w_j γ_j d^*(𝐱_i^t ,𝐱_c_k^j^t) = k'

holds almost surely for all i=1,…,N and all t≥max{T,T(N)}. 
Theorem <ref> is proved.

The next part involves discussion of the complexity costs of the above two algorithms. 


  * For offline setting, our Algorithm <ref> requires N(N-1)/2 calculations of d^*, against κ N calculations of d in the offline algorithm in <cit.>. In each d^*, the matrices distance ρ^* consists of m_n^2 calculations of Euclidean distances. Then iterating over m,l in d^* we see that at most 𝒪(nm_n^3) computations of Euclidean distances, against 𝒪(nm_n/|log s|) computations of d̂ for the offline algorithm in <cit.>, where

    s=min_ X_i^(1)≠ X_j^(2)
     i∈{1,…,n_1};j∈{1,…,n_2}|X_i^(1)-X_j^(2)|.

It is known that efficient searching algorithm can be utilized to determine s,  with at most 𝒪(nlog(n)) (n=min{n_1,n_2}) computations. Therefore our Algorithm <ref> is computationally competitive to the one in <cit.>.

  * For online setting, we can hold a similar discussion as in <cit.>, Section 5.1. There it shows the computational complexity of updates of d^* for both our Algorithm <ref> and the online algorithm in <cit.> is at most 𝒪(N(t)^2+N(t)log^3n(t)) (here we take m_n(t)=⌊log n(t)⌋). Therefore the overall difference of computational complexities between the 2 algorithms are reflected by the complexity of computing d^* and d (see Point 1).



 §.§ Efficient Dissimilarity Measure

<cit.> presented a set of three simple properties that a good clustering function should have: scale-invariance, richness and consistency. Further, he demonstrated that there is no clustering function that satisfies these properties at the meanwhile. He pointed out, as one particular example, that the centroid-based clustering basically does not satisfy the above consistency property (note that this is a different concept from our asymptotic consistency). In this section we show that, although the consistency property is not satisfied, there exists some other criterion of efficiency of dissimilarity measure in a particular setting. It is the so-called efficient dissimilarity measure.


Assume that the samples S={𝐱(ξ): ξ∈ℋ} (ℋ⊂ℝ^q for some q∈ℕ), meaning that all the paths 𝐱(ξ) are indexed by a set of real-valued parameters ξ. Then a clustering function is called efficient if its dissimilarity measure d satisfies that, there exists c>0 so that for any 𝐱(ξ_1),𝐱(ξ_2)∈ S,

    d(𝐱(ξ_1),𝐱(ξ_2))=cξ_1-ξ_2,

where · denotes some norm defined over ℝ^q.

Mathematically, efficient dissimilarity measure is a metric induced by some norm. Clustering processes based on efficient dissimilarity measure will then be equivalent to clustering under classical distances in ℝ^q, such as Euclidean distance, Manhattan distance, or Minkowski distance. The latter setting has well-known advantages in cluster analysis. For example,  Euclidean distance performs well when deployed to datasets that include compact or isolated clusters <cit.>; when the shape of clusters is hyper-rectangular <cit.>, Manhattan distance can be used; Minkowski distance, including Euclidean and Manhattan distances as its particular cases, can be utilized to solve clustering obstacles <cit.>. There is a rich literature on comparing the above three distances to each other through discussing of their advantages and inconveniences. We refer to <cit.> and the references therein.

In the next section we present an excellent example, to show how to improve the efficiency of our consistent algorithms, for clustering self-similar processes with wide-sense stationary ergodic increments.


§ SELF-SIMILAR PROCESSES AND LOGARITHMIC TRANSFORMATION


In this section we introduce a non-linear transformation of the covariance matrices in d^*, in order to improve the efficiency of clustering. This transformation is based on logarithmic function. We use one example to explain how this transformation works. We show this transformation maps d^* to some covariance-based dissimilarity measure similar to an efficient one, when applied to clustering self-similar processes.


A process X^(H)={X_t^(H)}_t∈ T (e.g., T=ℝ or ℤ) is self-similar with index H∈(0,1) if, for all n∈ℕ, all t_1,…,t_n∈ T,  and all c≠0 such that ct_i∈ T (i=1,…,n), 

    (X_t_1^(H),…,X_t_n^(H))law=(|c|^-HX_ct_1^(H),…,|c|^-HX_ct_n^(H)).


It can be shown that a self-similar process has necessarily zero mean and its  covariance structure is indexed by its self-similarity index H, in the following way <cit.>.


Let {X_t^(H)}_t∈ T be a zero-mean self-similar process with index H∈(0,1) and  with wide-sense stationary ergodic increments. Assume 𝔼|X_1^(H)|^2<+∞, then for any s,t∈ T,

    ℂov(X_s^(H), X_t^(H))=𝔼|X_1^(H)|^2/2(|s|^2H+|t|^2H-|s-t|^2H).


The corollary below follows.

Let {X_t^(H)}_t∈ T be a zero-mean self-similar process with index H and weakly stationary increments. Assume 𝔼|X_1^(H)|^2<+∞. For h>0 small enough, define the increment process Z_h^(H)(s)=X_s+h^(H)-X_s^(H), then for s,t∈ T such that s-t≥ h, we have

    ℂov(Z_h^(H)(s),Z_h^(H)(t) )=𝔼|X_1^(H)|^2/2((s-t-h)^2H+(s-t+h)^2H-2(s-t)^2H).


Applying three times the mean value theorem to (<ref>) leads to 

    ℂov(Z_h^(H)(s), Z_h^(H)(t) )=H𝔼 |X_1^(H)|^2((v_1^(H))^2H-1-(v_2^(H))^2H-1)h
          =H(2H-1)𝔼 |X_1^(H)|^2(v^(H))^2H-2h,

for some v_1^(H)∈(s-t,s-t+h), v_2^(H)∈(s-t-h,s-t) and v^(H)∈(v_2^(H),v_1^(H)).  We see that the item ℂov(Z_h^(H)(s),Z_h^(H)(t)) is a non-linear function of H. Next we would find a function g such that g(ℂov(Z_h^(H)(s),Z_h^(H)(t))) is linearly dependent of H. To this end we introduce the following log^*-transformation: for x∈ℝ, define

    log^*(x):=(x)log|x|={[   log(x)        ;; -log(-x)        ;;        0        . ].

Introduction to log^*-transformation is driven by the following 2 motivations:


  Motivation 1 The log^* function transforms the current dissimilarity measure to the one which linearly depends on its variable H.

  Motivation 2 The value log^*(x) preserves the sign of x, which leads to the consequence that larger distance between x,y yields larger distance between log^*(x) and log^*(y). 

Applying log^*-transformation to the covariances of Z_h^(H) given in (<ref>), we obtain

    log^*(ℂov(Z_h^(H)(s),Z_h^(H)(t)))
          =(2H-1)((2H-2)log v^(H)+log h+log(H|1-2H|𝕍ar(X_1^(H)))).

When v^(H) and h are small the items log v^(H) and log h are significantly large so log(H|1-2H|𝕍ar(X_1^(H))) becomes negligible. Thus we can write

    log^*(ℂov(Z_h^(H)(s),Z_h^(H)(t)))≈(2H-1)((2H-2)log v^(H)+log h).

In conclusion, 


  * When H_1,H_2∈(0,1/2] or H_1,H_2∈[1/2,1), the item log^*(ℂov(Z_h^(H)(s),Z_h^(H)(t))) is approximately linear on H∈(0,1/2] or on H∈[1/2,1). 

Using the approximation log v^(H_1)≈log v^(H_2) for H_1,H_2∈(0,1/2] or H_1,H_2∈[1/2,1), we have

    log^*(ℂov(Z_h^(H_1)(s),Z_h^(H_1)(t)))-log^*(ℂov(Z_h^(H_2)(s),Z_h^(H_2)(t)))
          ≈ 2(2H_1-1)(H_1-H_2)log v^(H_1).


  * When H_1∈(0,1/2] and H_2∈(1/2,1), log^*(ℂov(Z_h^(H)(s),Z_h^(H)(t))) turns out to be relatively large, because we have 

    log^*(ℂov(Z_h^(H_1)(s),Z_h^(H_1)(t)))-log^*(ℂov(Z_h^(H_2)(s),Z_h^(H_2)(t)))
          ≈-(2H_1-2)log v^(H_1)-(2H_2-2)log v^(H_2)
          ≥ 2(2-H_1-H_2)min{log v^(H_1),log v^(H_2)}.


Taking advantage of the above facts we define the new empirical covariance-based dissimilarity measure (based on the definition (<ref>)) to be

    d^**(𝐳_1,𝐳_2):=∑_m= 1^m_n∑_l= 1^n-m+1 w_m w_l ρ^* (ν^**(Z^(H_1)_l… n,m), ν^** (Z^(H_2)_l… n,m)),

where ν^**(Z^(H_1)_l… n,m) is the empirical covariance matrix of Z_h^(H_1),  ν^*(Z^(H_1)_l… n,m), with each of its coefficients transformed by log^*: let M={M_i,j}_i=1,…,m; j=1,…,n be an arbitrary real-valued matrix, define

    log^*M:={log^*M_ij}_i=1,…,m; j=1,…,n.

Then we have 

    ν^**(Z^(H_1)_l… n,m):=log^*(ν^*(Z^(H_1)_l… n,m)).

Now given 2 wide-sense stationary ergodic processes X^(1), X^(2), we choose {w_j}_j∈ℕ to satisfy

    ∑_m,l=1^∞  w_m w_l ρ^*(log^*(V_l,l+m-1(X^(1))),log^*(V_l,l+m-1(X^(2)))<+∞,

where we denote by

    V_l,l+m-1(X^(1)):=ℂov(X_l^(1),…,X_l+m-1^(1)).

Then define the log^*-transformation of the covariance-based dissimilarity measure to be

    d^**(X^(1),X^(2)):=∑_m,l=1^∞  w_m w_l ρ^*(log^*(V_l,l+m-1(X^(1))),log^*(V_l,l+m-1(X^(2))).

Using the fact that log^* is continuous over ℝ\{0} and the weak ergodicity of Z_h^(H), we have the following version of ergodicity:

    d^**(𝐳_1,𝐳_2)d^**(Z_h^(H_1),Z_h^(H_2)).

Unlike d^*, the dissimilarity measure d^** is approximately linear with respect to the self-similarity index H.  Indeed, it is easy to see that

    d^**(𝐳_1,𝐳_2)∼{[    |H_1-H_2|<1,               ;; 2(2-H_1-H_2)>1,               , ].

where H_1,H_2 correspond to the self-similarity indexes of X^(H_1),X^(H_2) respectively. In fact, from (<ref>) we can say that d^** satisfies Definition <ref> in the wide sense: it is approximately linearly dependent of |H_1-H_2| when H_1,H_2 are in the same group out of (0,1/2] and [1/2,1); it is approximately larger than |H_1-H_2| when H_1,H_2 are in different groups out of (0,1/2] and [1/2,1). This fact allows our asymptotically consistent algorithms to be more efficient when clustering self-similar processes with weakly stationary increments, having different values of H. In Section <ref> we provide an example of clustering using our consistent algorithms with and without the log^*-transformation, when the observed paths are from a well-known self-similar process with stationary increments – fractional Brownian motion.



§ SIMULATION AND EMPIRICAL STUDY
 
This section is devoted to applying our clustering algorithms to several synthetic data and real-world data. It is worth noting that, in our statistical setting, the auto-covariance functions are supposed to be unavailable, then the prior choice of the weights w_j presents some trade-off between the convergence of the dissimilarity measure and practical application. On one hand, low rate of convergence (e.g. w_j=1/j(j+1)) risks to a divergent dissimilarity measure d^* (see (<ref>)). On the other hand, high rate of convergence (e.g., w_j=1/j^3(j+1)^3) will only make use of some first observations in the sample paths. We believe that the first issue is a minor one in practice, because for most of the wide-sense stationary ergodic processes (especially Gaussian) taking w_j=1/j(j+1) can lead to convergent d^*. Also, in practice, instead of (<ref>) it is fine to regard

    d^*(X^(1),X^(2)):= ∑_m,l = 1^N w_m w_l ρ^*(V_l,l+m-1(X^(1)),V_l,l+m-1(X^(2))),

for some N large enough.

Therefore, through this entire section we take w_j=1/j(j+1) and m_n=⌊log n⌋ (recall that ⌊·⌋ denotes the floor number) in the covariance-based dissimilarity measure d^*. Next we explain how to prepare offline and online datasets in this simulation study.


  Offline dataset simulation: For each scenario, we simulate 5 groups of sample paths, each consists of 10 paths with length N(t)=5t, for the time steps  t=1,2,…,50. Algorithm <ref> is performed over 100 such scenarios, and the misclassification rate is calculated.

  Online dataset simulation: For each scenario, we simulate 5 groups of sample paths. Let the total number of sample paths be N(t) = 30 + ⌊ (t-1)/10 ⌋ at each time step t. That is, there are 6 sample paths in each of the 5 groups when t=1. And the number of sample paths in each group will increase by 1 once the time t increases by 10. For i=1,2,…, the ith sample path in each group has length n_i(t) = 5[t-(i-6)^+], where x^+ = max(x,0). 

We then apply the proposed clustering algorithms to both offline and online settings, and determine their corresponding misclassification rates. These misclassification rates are utilized to intuitively illustrate the asymptotic consistency of our clustering algorithms, or to compare the performances of our clustering approaches to other ones.  Recall that the misclassification rate (i.e. mean clustering error rate, see Section 6 in <cit.>) is obtained by dividing the number of misclassified paths by the total number of paths per scenario, then average all these fractions:

    p:=avg(# /# ).

More precisely, let (C_1,…,C_κ) denote the ground truth clusters of the N sample paths x_1,…, x_N. We define the ground truth cluster labels by

    L_k = (k,…,k)_# C_k  ,  k=1,…,κ.

Let (l_1,…,l_N) denote the cluster labels of (x_1,…,x_N) output by some clustering approach. Then the misclassification rate p of this approach is computed by

    p=min_σ∈ S_κ
    (π_1,…,π_N)=(L_σ(1),…,L_σ(κ))∑_i=1^N 1_{π_i l_i}/N,

where S_κ denotes the group of all possible permutations over the set {1,…,κ}.

For example, in one scenario of 7 sample paths, if the ground truth cluster labels of (x_1,…,x_7) satisfy

    (L_1, L_2, L_3)= ((1,1),(2),(3,3,3,3)),

while the clustering algorithm output cluster labels
 corresponding to (x_1,…,x_7) are given by 

    (l_1,…,l_7) = (2,1,1,2,3,2,1),

then according to Eq.  (<ref>), the misclassification rate is 4/7. This can be explained as, at least 4 changes of labels are needed to let the output cluster labels match that of the ground truth ones (1,1,3,2,2,2,2):

    l_1← 1; l_3← 3; l_5← 2; l_7← 2.

We provide the implementation of the misclassification rate (see Eq. (<ref>)) in MATLAB  publicly online as misclassify_rate.m [<https://github.com/researchcoding/clustering_stochastic_processes/blob/master/misclassify_rate.m>.].



 §.§ Clustering Non-Gaussian Discrete-time Stochastic Processes


In <cit.> a simulation study on a non-Gaussian  strictly stationary ergodic discrete-time stochastic process (see also <cit.>) has been performed. Since this process has finite covariance structure, it is also wide-sense stationary ergodic. As a result we can test our clustering algorithms over the same dataset and compare their performances to the ones in <cit.>. Recall that this process {X_t}_t∈ℕ is  generated in the following way. Fix some irrational-valued parameter α∈ (0,1).


  Step 1. Draw a uniform random number r_0 ∈ [0,1]. 

  Step 2. For each index i=1,2,…,N:


  Step 2.1. Define
r_i = r_i-1 + α - ⌊ r_i-1 + α⌋.

  Step 2.2. Define
X_i = 

1        when  r_i > 0.5, 

0        otherwise.


We simulate 5 groups of sample paths {X_t}_t∈ℕ indexed by the irrational values α_1 = 0.31..., α_2 = 0.33..., α_3 = 0.35..., α_4 = 0.37..., α_5 = 0.39... (α_i, i=1,…,5, each is simulated by a longdouble with a long
mantissa, see <cit.>), respectively. 


  §.§.§ Offline Dataset

We demonstrate the asymptotic consistency of Algorithm <ref> by conducting offline clustering on the simulated offline datasets of {X_i}_i∈ℕ. 

The valid blue line in Fig. <ref> illustrates the asymptotic consistency of Algorithm <ref> through the fact that its misclassification rate decreases as time t increases. Compared to the simulation study over the same dataset in <cit.>, the misclassification rate provided by our proposed algorithm converges at a comparable speed (see Figure 2 in <cit.>), even though Algorithm <ref> aims to cluster covariance structures but not process distributions. 

The dot-dashed red line in Fig. <ref> presents the performance of Algorithm <ref> and compares its misclassification rates with the ones from Algorithm <ref>. Applied to offline dataset, the offline algorithm's misclassification rates are consistently lower than the online algorithm, i.e., the offline dataset clustering algorithm performs better than the online dataset clustering algorithm, when dealing with offline datasets.






  §.§.§ Online Dataset

In our simulated online datasets the number of sample paths and the length of each sample path increase as t increases. This type of setting is mimicking the situation such as modeling financial asset prices, where new assets are launched at each time step. The offline and online clustering algorithms are applied at each time t with 100 runs, their misclassification rates at each time t are then obtained.
 
Fig. <ref> compares the misclassification rates of offline algorithm and online algorithm applied to the online dataset described above. The periodical pattern, that misclassification rate increases per 10 time steps using offline algorithm, matches the timing of adding new observations. That is, the misclassification rate spikes whenever new observations are obtained. We observe that the misclassification rate of the online algorithm is overall lower than that of offline algorithm in this dataset, reflecting the advantage of online algorithm against the offline one in the case where new observations are expected to occur. It is worth pointing out that our online setting is different from the one in <cit.>, therefore the two clustering results are not comparable.

Finally, all the codes in MATLAB that reproduce the main conclusions in this subsection can be found publicly online[<https://github.com/researchcoding/clustering_WSSP_with_cov_distance>.].




 §.§ Clustering Fractional Brownian Motions
 
In this section, we present the performance of proposed offline (Algorithm <ref>) and online (Algorithm <ref>) methods, on a synthetic dataset sampled from continuous-time Gaussian processes. The wide-sense stationary ergodic processes that we choose are the first order increment processes of fractional Brownian motions (see <cit.>). Denote by {B^H(t)}_t≥0 a fractional Brownian motion with Hurst index H∈(0,1). It is well-known that B^H is a zero-mean self-similar Gaussian process with self-similarity index H and with covariance function

    ℂov(B^H(s),B^H(t))=1/2(s^2H+t^2H-|s-t|^2H), .

Fix h>0, define its increment process (with time variation h) to be

    Z_h^(H)(t)=B^H(t+h)-B^H(t), .

Z_h^(H) is also called fractional Gaussian noise. Using the covariance function (<ref>) we obtain the auto-covariance function of Z_h^(H) below: for τ≥0,

    γ(τ)=ℂov(Z_h^(H)(s),Z_h^(H)(s+τ))=1/2(|τ+h|^2H+|τ-h|^2H-2|τ|^2H).

Recall that for stationary Gaussian processes such as Z_h^(H), the strict ergodicity can be fully expressed in the language of its auto-covariance function γ, i.e., the following result <cit.> provides a sufficient and necessary condition for a stationary Gaussian process to be strictly ergodic.


A continuous-time Gaussian stationary process X is strictly ergodic if and only if

    lim_t→∞1/t∫_0^t|γ_X(u)| u=0,

where γ_X denotes the auto-covariance function of X.

In view of (<ref>) we can deduce that the auto-covariance function γ of Z_h^(H) satisfies (<ref>). This together with Theorem <ref> yields that Z_h^(H) is second-order strict-sense stationary ergodic, so it is also wide-sense stationary ergodic.

To test our algorithms we simulate κ=5 groups of independent fractional Brownian paths, with the ith group containing 10 paths as {B^H_i(1/n),…,B^H_i((n-1)/n),B^H_i(1)}, for the self-similarity indexes

    H_1=0.3, H_2=0.4, …, H_5=0.7.

Remark that clustering a zero-mean fractional Brownian motion B^H is equivalent to clustering its increments Z_1/n^(H)(t)=B^H(t+1/n)-B^H(t). These total number of 50 observed paths of Z_1/n^(H)(t), each of length 150, compose an offline dataset and an online one.  The clustering algorithms are applied to the dataset at each time step t. 100 runs are made to compute the misclassification rates. we use offline (RESP. online) dataset clustering algorithm to cluster offline (RESP. online) dataset. The purpose is to compare the the algorithms with and without log^*-transformations.

Fig. <ref> presents the comparisons of 2 algorithms: one is using the dissimilarity measure d^*, the other one is using the dissimilarity measure d^**,  based on the behavior of misclassification rates as time increases. We conclude that, both algorithms with and without the log^*-transformations are asymptotically consistent. However in both offline and online settings, the covariance-based dissimilarity measure algorithms with log^*-transformation (dashed red lines) have 30% lower misclassification rates on average than that of algorithms without log^*-transformation (solid blue lines). This simulation study proves the necessity of utilizing log^*-transformed covariance-based dissimilarity measure when the underlying observations have nonlinear, especially power based, covariance-based dissimilarity measure, such as observations sampled from self-similar processes.

The codes in MATLAB used in this subsection are provided publicly online[<https://github.com/researchcoding/clustering_stochastic_processes>.].






 §.§ Clustering AR(1) Processes: Non Strict-sense Stationary Ergodic
 
To show that our algorithms can be applied to clustering non strict-sense stationary ergodic processes, we consider a simulation study on the non-Gaussian AR(1) process {Y(t)}_t defined in Example 2, Eq. (<ref>).  We then conduct the cluster analysis with κ=5, and specify the values of a in  Eq. (<ref>) as

    a_1 = -0.4,  a_2 = -0.15,  a_3 = 0.1,  a_4 = 0.35,  a_5 = 0.6.

We mimic the procedure in Section <ref> to generate the offline and online datasets of {X(t)}_t. Fig. <ref> illustrates the consistent converging property of offline algorithm and online algorithm under different dataset settings. 

All the codes in MATLAB that reproduce the main conclusions in this subsection can be found publicly online[<https://github.com/researchcoding/clustering_nonGaussian_processes>.].









 §.§ Application to the Real World: Clustering Global Equity Markets



  §.§.§ Data and Methodology

In this section we apply the clustering algorithms to real-world datasets. The application involves in dividing equity markets of major economic entities in the world into different subgroups. In financial economics, researchers usually cluster global equity markets according to either geographical region or the development stage of the underlying economic entities. The reasoning of these clustering methods is that entities with less geographical distance and closer development level involve in more bilateral economic activities. Impacted by similar economic factors, entities with less “distance” tend to have higher correlation in stock market performance. This correlation then measures the level of “comovement” of stock market indexes on global capital market. 

However, the globalization is breaking the barriers of region and development level. For instance, in 2016 China became the largest trader partner with the U.S. (besides EU)[Source: U.S. Department of Commerce, Census Bureau, Economic Indicators Division.]. China is not a regional neighbor of the U.S., and is categorized as a developing country by World Bank, in opposite to the U.S. as a developed country. 

We cluster the equity markets in the world according to the empirical covariance structure of their performance, using Algorithms <ref> and <ref> as purposed in this paper. Then we compare our clustering results with the traditional clustering methodologies. The index constituents of MSCI ACWI (All Country World Index) are selected as the sample data. Each of the observations is a sample path representing the historical monthly returns of underlying economic entities. Through empirical study it is proved that these indexes returns exhibit the long memory path feature hence they can be modeled by self-similar processes such as fractional Brownian motions (see e.g. <cit.>). Therefore similar to Section <ref> we may cluster the increments of the indexes returns with the log^*-transformed dissimilarity measure d^**.  MSCI ACWI is the leading global equity market index and has $3.2 billion in underlying market capitalization[As of June 30, 2017, as reported on September 30, 2017 by eVestment, Morningstar and Bloomberg.]. MSCI ACWI contains 23 developed markets, 24 emerging markets from 4 regions: Americas, EMEA (Europe, Middle East and Africa), Pacific and Asia. Table <ref> lists all markets included in this empirical study. We exclude Greece market due to its bankruptcy after the global financial crisis. 

We construct both offline and online datasets starting from different dates. For offline dataset we let it start from Jan. 30, 2009 to exclude the financial crisis period in 2007 and 2008. This is because, under global stock market crisis, the (downside) performance of equity market is contagious and thus blurs the cluster analysis. The online dataset starts on Jan. 31, 1989, which covers 1997 Asian financial crisis, 2003 dot-com bubble and 2007 subprime mortgage crisis. Another key feature is that 14 markets are added to the MSCI ACWI index (at different time) since 1989, including 1 developed market and 13 emerging markets. Therefore, the case where new time series are observed is handled in online dataset. 

[htbp]
The categories of major equity markets in the MSCI ACWI (All Country World Index). There are 23 markets from developed economic entities, and 24 markets from emerging countries or areas. The geographical clustering contains Americas, EMEA (Europe, Middle East and Africa), Pacific and Asia. 



 3cDeveloped Markets         3cEmerging Markets 

 (r)1-3 (r)5-7
Americas	   	Europe & Middle East	   	Pacific	       	Americas	   	Europe & Middle East & Africa	   	Asia	

 (r)1-3 (r)5-7
Canada	   	Austria	   	Australia    	   	Brazil	   	Czech Republic	   	China (Mainland)	

USA	   	Belgium	   	 Hong Kong	       	Chile	   	Greece	   	India	

	   	Denmark	    Japan	   	    Colombia	   	Hungary	   	Indonesia	

	   	Finland	   	New Zealand    	   	Mexico	   	Poland	   	Korea	

	   	France	   	Singapore    	   	Peru	   	Russia	   	Malaysia	

	   	Germany	   		       		   	Turkey	   	Pakistan	

	   	Ireland	   		       		   	Egypt	   	Philippines	

	   	Israel	   		       		   	South Africa	   	Taiwan	

	   	Italy	   		       		   	Qatar	   	Thailand	

	   	Netherlands	   		       		   	United Arab Emirates	   		

	   	Norway	   		       		   		   		

	   	Portugal	   		       		   		   		

	   	Spain	   		       		   		   		

	   	Sweden	   		       		   		   		

	   	Switzerland	   		       		   		   		

	   	United Kingdom	   		       		   		   		

 

 
Source: MSCI ACWI (All Country World Index) market allocation. <https://www.msci.com/acwi>.







  §.§.§ Clustering Results

We compare the clustering outcomes of both offline and online datasets with separations suggested by region (4 groups) and development level (2 groups). The factor with the lowest misclassification rate is proved to be the corresponding factor that contributes to increase covariance-based dissimilarity measure the most. In other words, this corresponding factor leads to the clustering of stock markets with the most significant impact. 

Table <ref> shows that the misclassification rates for development levels are significantly and consistently lower than that of geographical region, for both algorithms (offline and online algorithms) and datasets (offline and online datasets). The clustering results seem to infer that the geographical distance is less dominating than the development level of underlying economic entities, when analyzing different groups of equity markets.



The global minimum of the misclassification rate occurs when we use online algorithm on offline dataset. Table <ref> presents the detailed clustering outcome under this circumstance. In each group, the correctly and incorrectly categorized equity markets are listed respectively. For instance, China (Mainland) market is correctly categorized along with other emerging market. Meanwhile Austria market, though being developed market in MSCI ACWI, is categorized to the group where most of the equity markets are emerging markets. The misclassified markets in the emerging group are Austria, Finland, Italy, Norway and Spain markets. The misclassified markets in the developed group are Malaysia, Philippines, Taiwan, Chile and Mexico markets. These empirical results thus suggest that several capital markets have irregular post-crisis performance which blurs the barrier between emerging and developed markets.



The contribution of this real-world dataset cluster analysis is two-fold. First, we explored and determined the principal force that brings structural difference in global capital markets, which potentially predicts the comovement pattern of future index performance. Second, we provided new evidence on the impact of globalization on breaking geographical barriers between economic entities.




§ CONCLUSION AND FUTURE PERSPECTIVES

Inspired by <cit.>, we introduce the problem of clustering wide-sense stationary ergodic processes. A new covariance-based dissimilarity measure is proposed to obtain asymptotically consistent clustering algorithms for both offline and online settings. The recommended algorithms are competitive for at least two reasons: 


  * Our algorithms are applicable to clustering  a wide class of stochastic processes, including any strict-sense stationary ergodic processes whose covariance structures are finite. 

  * Our algorithms are efficient enough in terms of their computational complexity cost. In particular, a so-called log^*-transformation is introduced to improve the efficiency of clustering, for self-similar processes.  

The above advantages have been supported through the simulation study on  non-Gaussian discrete-time processes, fractional Brownian motions, non-Gaussian non strict-sense stationary ergodic AR(1) processes, and a real-world application: clustering global equity markets. The implementations in MATLAB of our clustering algorithms are provided publicly online.


Finally we note that, the clustering framework proposed in our paper focuses on the cases where the true number of clusters κ is known. The case for which κ is unknown is still open and left to future research. Another interesting problem is that, many stochastic processes are not wide-sense stationary but they get a tight relationship with the wide-sense stationarity. For example, a self-similar process does not necessarily have wide-sense stationary increments, but their Lamperti transformations are strict-sense stationary <cit.>; locally asymptotically self-similar processes are generally not self-similar but their tangent processes are self-similar <cit.>. Our cluster analysis sheds light on clustering the above processes. These topics can be left for future research. 
 






spbasic


