## Agglomerative Likelihood Clustering

Lionel Yelibi ∗ and Tim Gebbie †

Department of Statistical Sciences, University of Cape Town, Rondebosch 7701, South Africa (Dated: October 6, 2021)

We consider the problem of fast time-series data clustering. Building on previous work modeling the correlation-based Hamiltonian of spin variables we present an updated fast non-expensive Agglomerative Likelihood Clustering algorithm (ALC). The method replaces the optimized genetic algorithm based approach (f-SPC) with an agglomerative recursive merging framework inspired by previous work in Econophysics and Community Detection. The method is tested on noisy synthetic correlated time-series data-sets with built-in cluster structure to demonstrate that the algorithm produces meaningful non-trivial results. We apply it to time-series data-sets as large as 20,000 assets and we argue that ALC can reduce compute time costs and resource usage cost for large scale clustering for time-series applications while being serialized, and hence has no obvious parallelization requirement. The algorithm can be an effective choice for state-detection for online learning in a fast non-linear data environment because the algorithm requires no prior information about the number of clusters.

PACS numbers: 05.10.Ln, 75.10.Nr, 89.65.Gh

## I. INTRODUCTION

Summarizing large quantities of fast real-time feature time-series data that is sampled from environments with unknown dynamically evolving and non-linear interactions requires some sort of unsupervised learning 1 . The problem of unsupervised statistical learning in the context of financial market data was explored in prior work [47] where the ability to quickly approximate superparamagnetic cluster configurations [2, 11] from data was shown. Concretely, it was shown that the proposed algorithm recovered the correct super-paramagnetic cluster configurations near the entropy maxima. Previous cases studies include data clustering of stocks [47] and gene data in [12], temporal states of financial markets [23], and state-detection for adaptive machine learning in trading [15]. There is an endless variety of potential use-cases for this type of fast big-data clustering technology. The key problem with previous implementations related to the compute times for very large data-sets.

Building on prior work we propose and demonstrate an alternative to fast Super-Paramagnetic Clustering (fSPC) [47] using a modern and streamlined implementation of the 'Merging Algorithm' first suggested in [12], one that can recover the same or better cluster configurations but with significantly reduced compute times. The algorithm is in spirit similar to network science community detection algorithms [4, 36].

The analysis of structure of stock market correlation matrices constitute a large body of study in econophysics, quantitative finance, and network science. Representing

∗ lionel.yelibi@alumni.uct.ac.za

† tim.gebbie@uct.ac.za

1 For a review on the last 20 years on financial market correlation based data clustering see [25], more generally on data clustering see [17], and particularly on time-series clustering see [1, 22].

the relationships between financial assets in terms of networks has become an increasingly useful approach to explore, gain insights, and visualize the collective impact of various aggregate dynamical process on assets; be it the impact of a market crashes, or that of the iterative application of a portfolio optimization. This is convenient because graphs can be represented by adjacency matrices, and correlation matrices can then be mapped to network data before they are processed; here passed through a community detection algorithm.

Network data is usually represented by sparse matrices, this can operationally be seen as the main difference between data clustering, and community detection - it is the difference in representation. In order to map correlation matrices into adjacency matrices a threshold is often applied [30, 41, 48]. The threshold value above which correlated nodes are linked or unlinked is arbitrary. However, this can add a layer of complexity in the study of graphs of financial assets, while also being an additional source of noise that has little to do with the dynamic processes governing the relationships between assets, or objects more generally. Our approach differs from threshold based approaches in that it is a faithful data clustering model which works with dense correlation matrices - this can make it appropriate for complex systems in noisy environments but without the arbitrariness of threshold choice.

Here, we again exploit the Noh Ansatz [31] and the Maximum Likelihood Estimation (MLE) approach first introduced by Giada and Marsili [11, 12]. We call the new algorithm Agglomerative Likelihood Clustering (ALC); it has the benefit of being less computationally expensive than the Parallelized Genetic Algorithms (PGAs) previously implemented in [15, 16, 47]. The performance enhancement arises from being able serialize the algorithm into a brute-force search across cluster configurations so as to avoid unnecessary computational expense of genetic algorithm used in, e.g. the f-SPC implementation.

The paper is as follows: in Section (II) we introduce the Giada-Marsili model, in Section (III) we describe the optimization algorithm, in Section (IV) we consider clustering of synthetic generated time-series data, Section (V) provides additional tools for when the signal-to-noise ratio is relatively low, Section (VI) explores the solutions recovered for complex correlation matrices with nested block structure, Section (VII) discusses and compares our new method execution runtime to the previous algorithm and HDBSCAN,this is followed by the discussion and conclusion in Section (VIII) that highlights the performance of the algorithm we have introduced.

## II. GIADA-MARSILI LIKELIHOOD MODEL

The spectral analysis of stock market correlation matrices provides an intuition taking the form of the Noh ansatz [31]: that there exists a hierarchical structure in financial markets where individuals stocks are subcomponents of larger groups of assets, and each asset's individual returns is influenced by the collection of assets it belongs to. This can be expressed in the form of a simple generative model [31]:

$$x _ { i } = f _ { i } + \epsilon _ { i } \quad \quad \ \ ( 1 ) \quad d \varepsilon$$

where x i are the stock's features, f i the cluster-related influence, and glyph[epsilon1] i the node's specific effect.

In turn, this generative model can be iteratively nested within a more complex set of independent hierarchies as expanded out through the unexplained noise terms in the presence of top-down and bottom-up causal relationships [40, 43]. Here we restrict ourselves to a single collection that is broken up into groups with a simple noise terms for unexplained externalities.

We consider a group of N observations embedded in a space with dimensionality D as the features, every observation is assigned a spin value. One version of the ansatz models the observation features such that

$$x _ { i } = g _ { s _ { i } } \eta _ { s _ { i } } + \sqrt { 1 - g _ { s _ { i } } ^ { 2 } } \epsilon _ { i } \, \text{ \quad \ \ } ( 2 ) \, \text{ \ \ } \text{clu} \, \text{ }$$

where x i is one feature, g s i the intra-cluster coupling parameter 2 , η s i the cluster-related influence, and glyph[epsilon1] i the observation's specific effect and measurement error. A covariance analysis yields additional terms such as n s the size of cluster s , and c s the intra-cluster correlation 3 .

We explicitly mention that n s &lt; c s &lt; n 2 s must be enforced: the lower bound is required because g s is undefined for values of c s ≤ n s , and the upper bound requires

2 The thermal average 〈 g s 〉 can be used to reconstruct data-sets sharing identical statistical features of the original time-series by using Eqn. (2) [11]

3 Here n s = ∑ N i =1 δ s ,s i , c s = ∑ N i =1 ∑ N j =1 C ij δ s ,s i δ s j ,s , and g s = √ c s -n s n 2 s -n s [11, 15].

a strict inequality because Eqn. (4) is undefined when c s = n 2 s . We introduce a Dirac delta function 4 to model the probability of observing data in a configuration S close to criticality[26]:

$$\stackrel { \kappa } { \underset { \mathfrak { s } _ { i } } { \mathfrak { d } } } } ^ { \kappa } \stackrel { \kappa } { \underset { \mathfrak { s } _ { j } } { \mathfrak { d } } } & P = \prod _ { d = 1 } ^ { D } \prod _ { i = 1 } ^ { N } \left \langle \delta ( x _ { i } - ( g _ { s _ { i } } \eta _ { s _ { i } } + \sqrt { 1 - g _ { s _ { i } } ^ { 2 } } \epsilon _ { i } ) ) \right \rangle \,. \quad ( 3 )$$

This joint likelihood is the probability of a cluster configuration matching the observed data for every observations, and for every feature. The log-likelihood derived from P can be thought of the Hamiltonian of this Potts system [44]:

$$L _ { c } = \frac { 1 } { 2 } \sum _ { s \colon n _ { s } > 1 } \left [ \ln \frac { n _ { s } } { c _ { s } } + ( n _ { s } - 1 ) \ln \frac { n _ { s } ^ { 2 } - n _ { s } } { n _ { s } ^ { 2 } - c _ { s } } \right ] \,. \quad ( 4 )$$

The sum is computed for every feature, and represents the amount of structure present in the data. The value of L c is indirectly dependent on spins via the terms n s and c s .

There are advantages to this method relative to most industry standard alternatives: First, that L c is completely dependent on C ij , and the dimensionality of the data-set only plays a part in computing C ij . Second, it is adaptive: the number of clusters is not given as an input, unlike K-MEANS [38] or similar algorithms. Clustering configurations are randomly generated, and that which maximizes L c provides us with the number of clusters, and their compositions. The importance of this later feature should not be lost as it removes a layer of unnecessary hyper-parameter tuning. As such the algorithm competes with the likes of DBSCAN [7, 20] and HDBSCAN [6, 27].

## III. AGGLOMERATIVE LIKELIHOOD CLUSTERING

A traditional perspective often use when considering clustering problems is to try and demarcate methods into either those that implement either top-down, or bottomup algorithms. Using this perspective, top-down methods are then thought of as divisive and consist in starting with a single cluster as initial condition and splitting (or partitioning) the graph in additional clusters iteratively while minimizing the cost. On the other hand, bottomup methods initially start with each observations in their own clusters, and proceed to merge them iteratively [29].

The so called 'Louvain' algorithm [4] is agglomerative and can then be argued to implement the later bottomup approach to 'community detection' on networks. It

4 Let y i = x i -( g s i η s i + √ 1 -g 2 s i glyph[epsilon1] i ) , and δ y ( ) a Dirac delta function of y which is 1 when y = 0, and 0 otherwise.

is in spirit very similar to the Merging Algorithm (MR) developed by Marsili and Giada in [12].

Our previously implemented methods were based broadly on efficient implementations that streamlined PGA optimization frameworks, e.g. those proposed in [15, 16, 47] allow for all sorts of mutations. However, these approaches can be sensitive to initial conditions because at every step a new generation of individuals is mutated, evaluated, and then a group of the best candidates survives until the next algorithm's iteration - this can be arbitrarily path dependent. The f-SPC algorithm [47] secured additional computational advantages within this PGA framework by being able to exclude a computationally expensive mutation constraint within the genetic programming framework used in [15, 16].

However, genetic algorithm based approaches have disadvantages, the key ones are discussed in Table (I) and relate to the need for : i.) ambiguous stopping or convergence criteria, ii.) random mutations, and iii.) some sort of parallelization to reduce compute times.

- (I)1 Convergence Criteria : Assuming the existence of multiple local maxima it tries to navigate around these 'sub-optimal' solutions on its way to a global maximum. However, there is no certainty and it is just assumed that the algorithm stops once a criteria is met - the algorithm is explicitly stochastic in convergence.
- (I)2 Random Mutations : Because the algorithm applies random mutations, the population size, the number and diversity of mutations, and the number of generations all have an impact on the final result - this can introduce path dependence.
- (I)3 Parallelization : The algorithm requires evaluating the entire mutated population at every iteration. This has both a computational and memory cost as it requires loading the data (i.e. the correlation or similarity matrix) on every worker. The Likelihood evaluation itself is inexpensive but multi-processing adds CPU-overhead. This can be mitigated by using GPUs as in [15].

TABLE I: Disadvantages of the PGA algorithms for the likelihood L c .

## A. Greedy Merging

To build a fast general bottom-up merging algorithm we again start with all N spins in N cluster, but we iteratively merge clusters in a greedy fashion.

The Giada-Marsili Merging Algorithm (MR) implementation requires computing the change in likelihood ∆ L c : We consider three clusters C 1 , C 2 , and C 3 with C 3 = C 1 + C 2 where the addition operator '+' means clusters C 1 and C 2 are merged. Marsili and Giada define two cases for ∆ L c [12]:

$$\mathfrak { I } _ { 1 } \quad \text{ Case 1 } \colon \Delta L _ { c } = L _ { c } ( C _ { 3 } ) - \max [ L _ { c } ( C _ { 1 } ), L _ { c } ( C _ { 2 } ) ] \ ( 5 )$$

$$\text{$2\colon \Delta L_{c}=L_{c}(C_{3})-[L_{c}(C_{1})+L_{c}(C_{2})]$} \quad ( 6 )$$

In Case 1, as described by Eq. (5), C 3 would be a better cluster than any of C 1 and C 2 . Here we chose to use the more restrictive definition, Case 2, as defined by Eq. (6). The key insight is to realize that Case 2 requires that the new merged cluster must be better than the combination of the two individual sub-clusters. We can iteratively exploit this by building an algorithm that performs a comprehensive grid search over the space of clusters.

To implement this we can modify Eq. (4) by removing the sum and only compute the likelihood of individual clusters:

$$L _ { c } = \frac { 1 } { 2 } \left [ \ln \frac { n _ { s } } { c _ { s } } + ( n _ { s } - 1 ) \ln \frac { n _ { s } ^ { 2 } - n _ { s } } { n _ { s } ^ { 2 } - c _ { s } } \right ] \,. \quad ( 7 )$$

The objective at every iteration is to then maximize ∆ L c over every possible move. The implementation we adopt to generate the moves is inspired by innovations made in community detection methods, i.e. community detection [18] algorithms such as the 'Louvain algorithm' [4].

Using this type of community detection we return to a bottom-up agglomerative approach to quickly enumerate candidate configurations using a likelihood method to determine the efficacy of the particular configuration. The agglomerative likelihood clustering framework is a general framework suited for any useful choice of likelihood function.

Here we will make the choice of the modified GiadaMarsili likelihood function, and then using Case 2 from the Merging Algorithm (Eqn. (6)). There are other choices of likelihood that could be made e.g a multifactor cluster specification, or the likelihood function for K-means. The framework is general - but the choice made here is specific to the spin Hamiltonian framework.

## B. Iterating cluster configurations

The implementation requires keeping track of both the correlation matrix, and the cluster configurations at each iteration of the algorithm. In order to retain a flexible but dynamic representation of the cluster configurations we introduce the idea of the tracker array (see Table (II)). This object will store the list of objects represented in each cluster during algorithm updates . 5 The algorithm is initialized with singleton clusters.

Once the initialization step in Table (II) is performed, we can move onto the actual optimization steps required,

5 In python the numpy computed correlation matrix can be stored to a python dictionary so that it can be easily and efficiently modified. The tracker can be stored as a list variable.

these are described in Table (III). This requires a cluster update step, computing the likelihood changes, updating the record of configurations using that tracker lists, updating the correlation matrix and then checking the convergence criteria.

- (II)1 The Correlation Matrix: C is stored in a structure that can be easily modifiable whenever new entries are needed, additional clusters created, or when previously clustered objects are removed.
- (II)2 Singleton Initialization : All objects start in their own clusters.
- (II)3 Tracker List : We create a list which stores lists of objects: each list represents a cluster, and the labels inside the list are cluster members. The variable n s is the cluster size.

TABLE II: The ALC initialization uses a singleton configuration, a correlation matrix stored in an efficient manner that can be easily modified during the algorithm implementation, and a tracker list that records the configurations.

- (III)1 Clustering : Pick an object at random in the tracker's labels, cluster it with all other objects, and store the resulting ∆ L c values.
- (III)2 Maximize Likelihood Change : Find out the highest ∆ L c ; and then if it is bigger than 0 continue to the next step. Else the object cannot be clustered, is removed from the list and the process restarts.
- (III)3 Tracker Update : if there exists a positive ∆ L c , create a new label for the new cluster and its content is the union of the two clusters we merge.
- (III)4 Correlation Matrix Update : The correlation values are the sums of the correlations of the two clusters merged. The self-correlation is the sum of the intra-cluster correlations.
- (III)5 Iterative Convergence : This process is repeated until ∆ L c is non-negative.

TABLE III: The main routine of agglomerative likelihood clustering algorithm for time-series clustering consists of merging clusters to find a maximum increase in likelihood (See pseudo-code in Appendix (1))

## IV. SYNTHETIC DATA

## A. Noise and uncorrelated data

First, we start testing our algorithm on random uncorrelated data 'white noise'. We generate uncorrelated data from the student-t distribution with 3 degrees of freedom with different data-set sizes and lengths. in Figure (1c) we can observe the probability distribution of cluster sizes recovered by the algorithm. Uncorrelated data yields clusters despite the non-existence of clusters by design. Additionally in Fig. (1a) we show the number of 'noise' clusters with increasing data-set size: The linear relationship between number of clusters and data-set size is trivial, however when the cluster number is normalized with N we observe that it in fact grows slower than N . The curve in Fig. (1b) could serve as the upper bound on our method's solutions.

The problem of identifying clusters when they are in fact not present in the data can be explained by estimation noise inherent to sample covariance estimation. The ground truth is the identity matrix, but due to finite size effects the estimated off-diagonal values are nonzero. The clusters sizes appear fat tailed on a semi-log plot, and the mode of the distribution is 5.

The implications are multiple: Small clusters of size 5 and below could be statistically indistinguishable from noise, while larger clusters are less likely to be spurious. Furthermore, true uncorrelated data should yield solutions of singletons which is not absolutely the case because as previously mentioned the number of 'noise' doesn't grow as fast as N .

## B. Clustering Synthetic Student-t data

Complex systems are characterized by their multi-scale dynamics in time and space. Financial markets are a good example of this where stock prices are recorded and aggregated on time-scales ranging from milliseconds to daily prices, if not weeks, months and years. There can be vastly different if not incompatible data generating processes operating on different time-scales. At the same time, markets are intrinsically non-stationary, which means data going far back in the past becomes stale, and this has a clear impact on most investment models performance.

Here we consider correlated synthetic data and generate Student-t distributed data-sets of 500, 1500, and 3000 time-series of variables lengths 60, 250, and 500 using the one-factor model in Eq. (B1) (See Appendix (B)). Gaussian mixtures are readily used within the machine learning literature to test data clustering. Here we simply use the student-t distribution because stock returns distributions are known to be fat-tailed [10, 13]. This stylized fact is some what captured by simulating with the student-t distribution; we would like to investigate the potential effect this may have on our algorithm.

(a) Average Number of clusters per clustering solution as a function of with Data-set Size. (b) Normalized Number of Cluster respect to the data-set Size as of cluster size for uncorrelated clusterfunction of the data-set Size. (c) Semi-log plot of the distribution less data-sets.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00951v4-figure-dc71c429-559b-4a7d-a400-45daa35eb84e.jpg)

FIG. 1: Descriptive Analysis of data-clustering of uncorrelated clusterless time-series data.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00951v4-figure-e08ec8f0-6edd-46ea-bb2d-2726f90f4fc8.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00951v4-figure-b970d8ad-8d0a-4942-871e-635084cdd827.jpg)

(a) Simulated : Cluster derived correlated time-series cumulative returns for 500 simulated assets over 1000 days.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00951v4-figure-084758ca-6008-48a0-9ba0-2dae413a3074.jpg)

- (b) Ground Truth : The true correlation matrix MST

(c) Estimated : The estimated correlation matrix MST

FIG. 2: Using synthetic correlated time-series data created with Table (2). In subplot a.) 500 Normalized time-series cumulative daily returns colored by cluster. And respectively in b.) and c.) The Minimum Spanning Trees of the true cluster configuration, and the estimated correlation matrix from the synthetic data. The colours represent the 10 clusters present in the data.

The data-sets of sizes: 500, 1500, and 3000 seem arbitrarily chosen. However, the size of the S&amp;P500 is roughly 500, and Indices of the Russell family: Russell 1000, 2000, and 3000 track up to 90% of the US stock market total market capitalization. Large portfolio of thousands of stocks exist for which equivalently large correlation matrices must be estimated. Correlation matrices are readily directly used in Markowitz's portfolio optimization [19, 49] or indirectly via clustering using methods such as hierarchical risk parity [35].

The intra-cluster coupling strengths g s , were arbitrarily selected to illustrate two limiting cases: 0.05 and 1,

TABLE IV: The Adjusted Rand Index (ARI) for the ALC and HDBSCAN cluster solutions for data-sets containing 10 clusters of correlated time-series. From left to right the size varies from 500 to 3000 while from top to bottom time-series length goes from 20 to 250. g s , the average intra-cluster correlation, is also changed from 0.05 to 1.

|              | ALC         | ALC   | ALC      | HDBSCAN   | HDBSCAN   | HDBSCAN   |
|--------------|-------------|-------|----------|-----------|-----------|-----------|
|              | N = 500 N = | 1500  | N = 3000 | N = 500   | N = 1500  | N = 3000  |
| g s = 0 . 05 |             |       |          |           |           |           |
| D = 20       | 0.04        | 0.06  | 0.03     | 0.03      | 0.04      | 0.00      |
| D = 60       | 0.07        | 0.07  | 0.07     | 0.03      | 0.03      | 0.02      |
| D = 250      | 0.19        | 0.14  | 0.14     | 0.09      | 0.08      | 0.09      |
| g s = 0 . 1  |             |       |          |           |           |           |
| D = 20       | 0.09        | 0.08  | 0.11     | 0.04      | 0.01      | 0.04      |
| D = 60       | 0.19        | 0.15  | 0.14     | 0.12      | 0.07      | 0.03      |
| D = 250      | 0.47        | 0.39  | 0.32     | 0.23      | 0.18      | 0.12      |
| g s = 0 . 3  |             |       |          |           |           |           |
| D = 20       | 0.27        | 0.24  | 0.24     | 0.22      | 0.17      | 0.12      |
| D = 60       | 0.58        | 0.49  | 0.41     | 0.40      | 0.31      | 0.20      |
| D = 250      | 0.90        | 0.84  | 0.76     | 0.66      | 0.60      | 0.35      |
| g s = 1      |             |       |          |           |           |           |
| D = 20       | 0.61        | 0.49  | 0.44     | 0.56      | 0.51      | 0.48      |
| D = 60       | 0.90        | 0.87  | 0.81     | 0.83      | 0.80      | 0.73      |
| D = 250      | 0.99        | 0.98  | 0.96     | 0.98      | 0.94      | 0.92      |

and two intermediate cases: 0.1, and 0.3. Here the value 0.3 is an approximation to the average stock market average correlation [34]. This is not the average intra-cluster correlation, which is what g s is; nevertheless, it gives a sense of how the algorithm would perform given real stock market data. The limiting values 0.05 and 1 serve to show how performance degraded ( g s = 0 05) when dealing with . low or weakly bound clusters, as opposed to perfectly correlated ones ( g s = 1). Given two time-series belong to different clusters, their pairwise correlation should be 0. However, this is not always the case because there is a noise-floor, and a lower g s will bring the clusters close to the noise floor. We have empirically confirmed in Table IV that this will deteriorate the algorithms performance. Here each of the size and length pairs was simulated a hundred times, and the Adjusted Rand Index (ARI) [37] between the ground truth and the model's output was averaged (See Table (IV)).

The time-series length were chosen with the noise floor effects in mind. Short time-series are expected to induce a higher noise floor, and spurious correlations. While longer time-series ameliorate this even though finite size effects remain. Towards this end we select the lengths: 20, 60, and 250. These are are approximations for the number of trading days in 1, 3, and 12 months of real world data. These lengths are adequate for data clustering on daily prices. It also follows that within a high frequency regime, one where measurements are abundant, correlation matrices are easier to estimate, and the algo- rithm should perform accordingly. However, in the highfrequency domain one is more directly concerned about the applicability of sample estimated correlations that do not deal with the impact of discretization [33].

## V. NOISE: RE-SAMPLING CLUSTERS AND BOOTSTRAP

In Sec. (IV B) we show that the estimation of correlation matrices induces statistical noise which significantly affects the performance of clustering algorithms. Here we provide a way to mediate the influence of noise using bootstrapping. The intuition arises from the relationship between the algorithm performance with respect to the signal to noise ratio Q = D N . The analysis of correlation matrices using Random Matrix Theory shows that they can be exactly estimated in the limit where N → ∞ , D →∞ and Q ≥ 1 and is maintained.

Despite the illusion of living in the so-called era of 'Big Data', in Finance, and more specifically in the context of financial markets, observables remain non-stationary and subject to extreme events, shocks and regime changes, all in the presence of strategic purposeful agents. Data quickly becomes obsolete. This means that the length of time-series used to estimate the correlation matrix C ij can be pragmatically shortened, or adaptively estimated, to better capture more relevant, time or regime dependent recent dynamics, but at the cost of increasing esti-

mation noise. As shown in Table (IV) this can have a severe impact on clustering solutions.

## A. Filtered correlation matrix

Here we consider cases in which the number of stocks N is large, but the number of realisations D is small, leading to Q ≤ 1 to violate the quality conditions. We then propose that a definitive filtered cluster membership matrix be constructed by sampling n in N stocks such that the quality condition on the n -sampled sub-problem is satisfied: q = D n ≥ Q . The routine that implements this is made explicit in Table (V). This introduces the idea of a filtered correlation matrix c ij , which is inspired from the original Potts 'spin-spin correlation function' [44].

Consider M random draws of n objects from a collection { 1 , ..., N } . First, from each sample we computed the pair-wise frequency of objects being drawn together, f ij :

$$f _ { i j } = \sum _ { m \in \mathcal { M } } \sum _ { i < j } \delta _ { i j }. \quad \quad ( 8 ) \quad \text{$\frac{1}{2}$ e\frac{1}{2}$} _ { t e r }$$

Second, we find the pair-wise frequency of objects begin clustered together, d ij :

$$d _ { i j } = \sum _ { m \in \mathcal { M } } \sum _ { i < j } \delta _ { s _ { i } s _ { j } }. \quad \quad ( 9 ) \quad \quad,$$

These can then be combined into a pair-wise normalised probability of clustered objects: p ij :

$$p _ { i j } = \frac { d _ { i j } } { f _ { i j } }. \quad \quad \quad ( 1 0 ) \quad \text{ as}$$

This is a measure of the probability that objects will being clustered together. This is then mapped to an ordinal codependency function on [0 1] using a threshold , ω :

$$p _ { i j } ^ { \omega } = \begin{cases} 1, & p _ { i j } - \omega > 0, \\ 0, & p _ { i j } - \omega \leq 0. \end{cases} \quad \begin{pmatrix} \ u e \\ ( 1 1 ) & \ u n e \\ \text{ind} \\ \text{nmr} \end{cases}$$

The final step is similar to the thresholding method used in the original Super-Paramagnetic Clustering (SPC) algorithm proposed by Blatt et al. [3] 6 . We note that any activation function could be used, and it is reminiscent of a logistics regression.

6 The original SPC uses simulated annealing, which via a repeated pairwise linkage validation process estimates a quantity similar to p ij (See. Eq. (10)) called the 'spin-spin correlation'

## (V)1 Hyper-parameters ( q n ω m , , , ):

Target signal-to-noise ratio: q

Sample size n computed from q

Correlation threshold: ω

Number of iterations: m

. . . .

(V)2 Resample Clusters: k = 1 to m

- 1. Compute N × N matrix f ( k ) ij ,
- 2. Cluster sample of n indices,
- 3. Compute N × N matrix d ( k ) ij .

(V)3 Probability matrix : Compute the N × N matrix p ij (Eqn. (10)).

(V)4 Ordinal Filtered Dependency matrix : Compute p ω ij (Eqn. (11))

(V)5 Final Configuration : The final cluster configuration can then be built from the ordinal matrix p ω ij acting as the adjacency matrix of a graph.

TABLE V: The bootstrapping routine for Agglomerative Likelihood Clustering (See Python script 'cluster resampling.py' [45]). The goal is to repeatedly cluster subsets of the data-set aggregated in d ij (Eq. (9)), transformed in normalized probability p ij (Eq. (10)) and filtered into ordinal filtered dependency p ω ij (Eq. (11)) from which the final cluster is extracted.

The n indices sampled are included in N , this means that the routine will cluster small subsets of the dataset for which the correlation matrix has a higher quality factor q than Q . The operation is repeated for a preset number of iterations m , and we use both matrices A and F to create the C ∗ . The threshold function is then applied to C ∗ to create a filtered adjacency matrix with disconnected components from is extracted the cluster solution: C F ∗ .

The number of maximum iterations m is arbitrarily chosen, in experiments we show that the convergence of the algorithm depends on n and indirectly on N . This means that larger data-sets require more sampling of the index space. We arbitrarily set the maximum iteration number: m = 2200.

Finally, we explored two threshold values: ω = 0 5 and . ω = 0 75 and we compare this to using the niave corre-. lation matrix. Although arbitrary, the intuition comes from [3] where the threshold used to cluster objects in the same group is a half: 0 5. .

The bootstrap method was tested on several correlated data-sets whose sizes and time-series length ranged from N = 1000 and D = 10 to N = 3000 and D = 30. The n sample size was chosen such that q = 0 1 given . Q = 0 01. . The output cluster solutions were compared against the ground truth using the ARI. The ARI is also used as a stopping condition: if ARI ≥ 0 9 the algorithm stops. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00951v4-figure-2cdff9b6-b0b8-47e7-aaf4-2d2cddabf0c4.jpg)

(a)

ARI for Bootstrap ALC using the threshold

ω

= 0 5.

.

(b)

ARI for Bootstrap ALC using the threshold

ω

= 0 75.

.

FIG. 3: The ARI is computed for the cluster solutions recovered using the bootstrap method and the ground truth. Data-sets size varies from 1000 to 3000, and the Signal-to-Noise ratio is fixed at q = 0 1 .

Figures (3a) and (3b) show thresholding with 0.5 is significantly worse than 0.75. In Fig. (3a) the solutions have yet to converge after 2000 iterations. Although they seem to be following an upward trend all trajectories are absolutely inferior to their Fig. (3b) counterparts after 500 iterations.

In Fig. (3b) we show the superior case. Given all data-sets have a fixed q = 0 1 we show that the bootstrap . method performance is correlated to the size of n with the larger sample sizes achieving better outcomes with less iterations. While the growth in performance is positively steep for the first 750 iterations it tends to slow down past that point for the 3 cases ( n = 100 150 200) after , , which the algorithm was stopped for reaching maximum iterations.

It is misleading to conclude that larger n is beneficial, this is because q was fixed while n was increasing. The issue relates back to Table (IV) where in the case of shortest time-series (i.e length 20) with data-sets size going from N =500, 1500, 3000 the respective ARI are 0.68, 0.54, and 0.48. The ARI is negatively correlated with the size of the data-set when time-series length is small. The bootstrap method applied to a similar case with N =2000, n = 200, and series length of 20 (See green curve in Fig. (3b)) yields an ARI higher than 0.8 after 750 iterations thus demonstrating the usefulness of the method when dealing with large data-sets and short time-series.

## VI. HIERARCHICAL BLOCK CORRELATIONS

ALC is efficient at clustering systems of correlated time-series generated using a single factor model when compared to HDBSCAN, but the likelihood model used is one explicitly derived using a single factor. Realistic correlation structures would include more complex timeseries models that incorporate multiple factors. In particular, an important related problem in financial timeseries clustering is that of clustering correlation matrices with a block hierarchical structure [24, 40, 41]. Here each block corresponds to a correlation cluster that should be recoverable along with its nested structure. Neither ALC nor HDBSCAN recover the dendrogram, and hence cannot be expected to recover the full nested structure without iterating recursively over each identified cluster until only singletons remain and retaining the configuration paths for these hierarchical iterations. For this reason it is useful to compare the difference between ALC and HDBSCAN on a block hierarchical problem as this will demonstrate more clearly how ALC deals with subclusters when they may exists within clusters.

Consider N stocks whose price return x i are influenced by a set of c hierarchical factors η , and glyph[epsilon1] the stocks individual effects; here η and glyph[epsilon1] are IID random variables ∼ N (0 , 1). The correlation matrix can be organized in nested clusters: h 1 ⊃ ·· · ⊃ h k ⊃ ·· · ⊃ h c , where h , 1 · · · , h c represent hierarchies at which factors η , 1 · · · , η c are sequentially added to the signal so that along the hierarchical path all stocks are correlated while lower hierarchies inherit all factors from upper ones. The picture is clearer when we consider a model of price return process x i :

$$h _ { 1 } \colon x _ { i } = \beta _ { 1 } \eta _ { 1 } + \alpha _ { 1 } \epsilon _ { i },$$

$$h _ { 1 } \colon x _ { i } = \beta _ { 1 } \eta _ { 1 } + \alpha _ { 1 } \epsilon _ { i }, & & ( 1 2 ) \\ & \vdots$$

$$h _ { k } \, \colon x _ { i } = \beta _ { 1 } \eta _ { 1 } + \cdots + \beta _ { k } \eta _ { k } + \alpha _ { k } \epsilon _ { i }, \text{ \quad \ \ } ( 1 3 )$$

.

.

.

$$h _ { c } \, \colon x _ { i } = \beta _ { 1 } \eta _ { 1 } + \cdots + \beta _ { k } \eta _ { k } + \cdots + \beta _ { c } \eta _ { c } + \alpha _ { c } \epsilon _ { i }. \, ( 1 4 )$$

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00951v4-figure-a955e8f8-c191-4fe9-b80b-35e788c71458.jpg)

(a) Correlation Matrix of a cluster with 5 hierarchies generated using Eqns. (12)-(14)

(b) Correlation Matrix of system of 3 clusters with 5 hierarchies generated using Eqns. (12)-(14)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00951v4-figure-a71299e4-042a-42cb-95b7-c9bdd8a0ca00.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00951v4-figure-a86472f0-0ae8-4a76-90e9-2fa7e31df833.jpg)

(c) Dendrogram of a cluster with 5 hierarchies generated using Eqns. (12)-(14) (d) Dendrogram of system of 3 clusters with 5 hierarchies generated using Eqns. (12)-(14)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00951v4-figure-8ac5247f-982f-4cd3-828b-4fe77ed4b5a8.jpg)

FIG. 4: Correlation matrices with nested clusters and their respective dendrograms. On the left is the 5-hierarchy model, and on the right 3 clusters with 5 nested hierarchies.

Here α k = 1 -∑ k j β j where k is the hierarchy index, α modulates the signal to noise ratio and can be arbitrarily set to any value by enforcing a constraint on the sum ∑ j β j to be equal to any value between 0 and 1. The β s are also chosen such that β j &lt; β j +1 this ensures that stocks belonging to hierarchy h k are more correlated to η k . There are many ways of parameterizing such nested linear models, and in fact this is just a special case of the more general model of Tumminello et al. [40].

system with 3 clusters estimated on timeseries of 250 observations ( i.e. 1 trading year) for α = 0 4 (as shown in . Fig. (4b)). Unlike our toy model, real correlation matrices are unlikely to be statically ordered in a way which can reveal their block structure which would trivialize clustering.

Using this model we demonstrate clustering on two hiearchical block correlation examples: a.) a correlation matrix of a cluster of 250 stocks with 5 hierarchies (50 stocks per hierarchy) in Fig. (4a), and b.) that of a

The dendrograms in Figures (4c) and (4d) created using the single linkage algorithm reveals the hierarchical structure present in the the correlation matrix in Figures (4a) and (4b). At first sight looking at the ordered correlation matrices one can distinguish the 5 hierarchies in Figures (4c) and (4a). Whereas in system containing multiples of these sets of nested clusters the apparent

structure is that of 3 clusters with potential sub-clusters.

The key question here is whether or not the solutions resulting from ALC and HDBSCAN contain the subclusters or the large clusters; i.e. do the results from the data analysis in Fig. (4a) contain 1 or 5 clusters, and those of Fig. (4b), 3 or 15 clusters? To answer this we have clustered these two datasets for several values of α to investigate how low correlation, and indirectly low cluster density, affect the algorithms' outputs. We show that for low α , the timeseries are mostly influenced by the factors, and HDBSCAN and ALC both recover 5 and 15 clusters respectively for the data in Fig. (4a) and (4b). The effect of noise is stronger for HDBSCAN, where the algorithm starts adding observations to the 'noise' category for datasets generated with α &gt; 0 2 (See Fig. . (6)) for data in Fig. (4a). While the ALC generated solutions start degrading from α &gt; 0 4. . Similarly for the data in Fig. (4b) HDBSCAN starts merging the nested clusters into 3 large clusters for α &gt; 0 2, . and for α &gt; 0 5 . the solution degrades even further (See Fig. (6)). Again, the ALC solutions are more resilient and start degrading for α &gt; 0 4 where nested clusters are merged and the 3 . cluster solution is never recovered.

This marks a distinct difference between ALC and HDBSCAN: ALC seems to prioritize smaller but more densely correlated clusters, whereas HDBSCAN default behavior is to merge clusters even if the resulting cluster is less dense (i.e. noisier). ALC seems to also be 20% more robust to noise. Finally for α &gt; 0 7 ALC solutions . are complete disordered partitions.

## VII. PERFORMANCE IMPROVEMENT

The appeal of ALC class of clustering algorithms is its ability to consistently outperform the f-SPC algorithm[47] while being generally robust to noise. Here we clustered data-sets of time-series using the synthetic model described in Appendix B with increasing size (N=50,100,200,300) and 10 clusters. Fig. 7a shows that ALC recovers clustering solutions of better quality if we use the likelihood as a quality function as was argued in [12]. Furthermore the algorithm runtime is roughly quadratic. This was estimated on synthetic data-sets 7 O N ( 1 97 . ) and on real data-sets 8 O N ( 2 11 . ) thus making it competitive for online learning problems. We compared the algorithm's runtime to the previous f-SPC and HDBSCAN [27]. Not only f-SPC's solutions have lower likelihood as previously mentioned but require much more computing power to achieve convergence in a significantly

7 Gaussian mixtures of 10 clusters of sizes ranging from 100 to 10000

8 We obtained CRSP US Mutual Funds Net Asset Values (NAV) sampled daily from 1998 to 2020 from which were sampled datasets of size from 100 to 10000

longer time . 9

Figure 7b demonstrates that future optimization algorithms cannot be expected to more expensive than ALC if they provide low quality solutions. Finally, HDBSCAN which we believe the natural alternative to ALC as a hierarchical and density based clustering method, is more optimized than our algorithm, and would thus be currently preferable for extreme large data-sets. However, there is some evidence that suggests that ALC can be considered as a serious alternative to HDBSCAN because under certain circumstances, when data is less dense, less correlated, and noisier, it outperforms in terms of ARI (See Table IV). The key point is that ALC looks to be competative both on quality and performance prior to any low level attempts at refining algorithm for optimised runtime performance.

## VIII. DISCUSSION AND CONCLUSION

We have presented an agglomerative algorithm capable of performing the maximization of the Giada-Marsili L c likelihood (See Eqn (4)). In prior work we built and proved a mechanism which maximizes the likelihood locally using Eqn. (7) instead of Eqn. (4) [47]. Here, instead of random moves, we perform a comprehensive search over all the possible combinations and select the optimal move at every iteration.

The algorithm is considerably faster than Markov Chain Monte Carlo, and Genetic Algorithm based solutions and achieves equal or better maxima. It only requires a correlation matrix as input, and its output is the clustering configuration which reflects the optimal number of clusters of correlated samples. It requires no prior information on the number of clusters. This feature may make the algorithm appropriate for state-detection for online learning in a fast-big-data environment.

We also present a data generation method of simulated correlated time-series based on the Noh ansatz given in Eqn. (2). Individual time-series are identified by their daily random effect glyph[epsilon1] while they remain subject to their cluster effect η , and the strength of that coupling between individual objects and their cluster is expressed by g s . This allowed us to modulate cluster densities, and study its impact on clustering quality and benchmark performance. We showed that even at low intra-cluster coupling values ( g s ≈ 0 05) the performance is poor, but . as g s increases and clusters densify the algorithm gets better faster than its competing alternative HDBSCAN.

9 For ALC and HDBSCAN, we are able to cluster larger datasets and extend the simulations up to N=10,000. All scripts were executed on a mobile Intel i7-CPU with 4 CPUs. This is stall generation dependent, and for a higher number it should be possible to obtain better results. This illustrates that f-SPC requires non-trivial parameterization, while ALC does not.

(a) The number of clusters recovered using HDBSCAN and ALC on data in Fig. (4a) as a function cluster specific and ALC on data in Fig. noise. ALC is more robust to cluster specific noise in the cific noise. hierarchical correlation model.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00951v4-figure-381f8d14-084a-44d7-8a7e-2966c964c3d9.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00951v4-figure-ac3d0a36-5359-4116-8d00-869112926d0f.jpg)

(b) The number of clusters recovered using HDBSCAN (4b) as function of cluster speALC is more robust to cluster specific noise in the hierarchical correlation model.

FIG. 5: The number of clusters as a function of the cluster specific noise level parameter α . As α increases the cluster specific noise increases. ALC degrades from α = 0 4 and at high levels of noise generates many singletons. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00951v4-figure-c733ce06-9f42-4748-9991-6a5023286406.jpg)

perhaps more efficient sampling schemes for the bootstrap step can mitigate this added cost [21]. The current Bootstrap method introduced in (V) uniformly samples a space of n in N variables which for large N becomes very expensive. Additionally, we demonstrated its use on synthetic data with a known ground truth. On real data, there is a clear need to define a stopping criteria for potential users (i.e. Convergence of the Likelihood L c .).

FIG. 6: The percentage of objects misclassified as noise by HDBSCAN applied to data in 'Cluster' Fig. (4a) and 'System' Fig. (4b)

When the signal-to-noise ratio is low it is possible under certain circumstances to mitigate the impact of correlation matrices noise by making use of bootstrapping in conjunction with ALC. This allows for better overall performance, but at the cost of significantly increasing the algorithm runtime. The increased computation related to bootstrap can be carried out in a batch format, or in parallel to the cluster algorithm itself. However, an ability to perform accurate clustering for relatively short timeseries (or equivalently lower dimensional data-sets) is valuable. We note that the exploration of alternate, and

When we design and cluster datasets using hierarchical factor models, both ALC and HDSCAN recover the underlying nested clusters. ALC is 20% more robust to noise for the synthetic datasets generated. In fact for large systems of such hierarchical clusters HDBSCAN has a tendency to merge clusters into giant ones. There is no single solution to this kind of clustering problem and it is left to the practitioner to determine whether clustering solutions with large but less dense clusters are more useful than solutions with small but denser clusters.

Potential further research could be achieved by using this algorithm in an online learning environment. We suspect in the case of financial markets that it may be possible to perform temporal clustering which would allow the analysis of the existing dynamics in financial markets states with increased reliability. Specifically for dynamic cluster analysis of financial markets around recorded extreme events 10 .

The current version of the algorithm is able to comfortably process data-sets of size up to 10,000 samples 11

10 Forthcoming work

11 Here tested on a Intel(R) Core(TM) i7-6700HQ CPU @ 2.60GHz

(a) Comparing cluster quality using the likelihood L c (See the models in Section II) of solutions obtained as the seconds data-set size is increased on a log-scale using f-SPC [47] SPC and ALC (Section III). ALC solutions are systematically are roughly quadratic, f-SPC significantly more expensive, of higher likelihood than ones obtained using f-SPC. (b) Comparing the algorithm runtime performance in as we increase data-set size on a log-scale for f[47], ALC (Section III), and HDBSCAN [27]. All while HDBSCAN outperforms ALC.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00951v4-figure-e3790ead-0e1b-4cf9-837b-f627cd37bf4d.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00951v4-figure-36ea8a6c-1934-41e8-a32f-9835e48d379e.jpg)

in less than 5 mins - enough time for a 'cup of coffee'. Crucially, it should also be noted that the Louvain algorithm exists in multiple faster versions than the original implementation use in the current algorithm [9, 32, 39] It feasible that one should be able to find further optimizations which will make it possible to cluster massively large dense correlation matrices similarly to what is currently done in the field of Network Science.

vates alternative modeling with distributions which capture additional stylized fact [10, 13].

The impact on distributional assumptions is also an important area of future work, the Giada-Marsili Likelihood of Section (II) was derived assuming IID Gaussian random variables. Using the same Noz Ansatz it should be possible to arrive at relate model extensions that assume different distributions . e.g. the Student-t distribution. Concretely, the Gaussian assumptions serves as a baseline for this type of model, but it is known that stock market log returns are at least fat tailed distributed (beside additional stylised facts such as long-memory, volatility clustering and the leverage effect). This moti-

[1] Aghabozorgi, S., Shirkhorshidi, A. S., and Wah, T. Y. (2015). Time-series clustering - A decade review. Information Systems , 53:16-38.

[2] Blatt, M., Wiseman, S., and Domany, E. (1996). Superparamagnetic Clustering of Data. Phys. Rev. Lett. , 76(18):3251-3254.

[3] Blatt, M., Wiseman, S., and Domany, E. (1997). Data Clustering Using a Model Granular Magnet. Neural Computation , 9(8):1805-1842.

Furthermore, different measures of similarity could be used. The covariance between two random variables is the standard similarity metric used to perform data clustering of financial market assets return series. Alternatives such as the information theoretic Mutual Information have previously been used directly on financial data time-series [8, 14].

Finally we believe it is possible to upgrade ALC with more efficient coding, and superior modifications to the optimization scheme itself similar to recent Louvain's implementations.

## IX. ACKNOWLEDGEMENTS

The authors thank Nic Murphy, Daniele Marinazzo, Gautier Marti and Unarine Singo for discussions and comments.

[4] Blondel, V. D., Guillaume, J.-L., Lambiotte, R., and Lefebvre, E. (2008). Fast unfolding of communities in large networks. Journal of Statistical Mechanics: Theory and Experiment , 2008(10):P10008.

[5] Bun, J., Bouchaud, J.-P., and Potters, M. (2017). Cleaning large correlation matrices: Tools from Random Matrix Theory. Physics Reports , 666:1-109.

[6] Campello, R. J. G. B., Moulavi, D., and Sander, J. (2013). Density-Based Clustering Based on Hierarchical Density Estimates. In Pei, J., Tseng, V. S., Cao, L., Motoda, H., and Xu, G., editors, Advances in Knowledge Discovery and Data Mining , Lecture Notes in Computer Science, pages 160-172, Berlin, Heidelberg. Springer.

- [7] Ester, M., Kriegel, H.-P., Sander, J., and Xu, X. (1996). A Density-based Algorithm for Discovering Clusters a Density-based Algorithm for Discovering Clusters in Large Spatial Databases with Noise. In Proceedings of the Second International Conference on Knowledge Discovery and Data Mining , KDD'96, pages 226-231. AAAI Press.
- [8] Fiedor, P. (2014). Networks in financial markets based on the mutual information rate. Phys. Rev. E , 89(5):052801.
- [9] Fontolan, F. (2020). Modularity Based Community Detection on the GPU . Bachelors Thesis, Universit` a Ca' Foscari Venezia.
- [10] Gabaix, X., Gopikrishnan, P., Plerou, V., and Stanley, H. E. (2003). A theory of power-law distributions in financial market fluctuations. Nature , 423(6937):267-270.
- [11] Giada, L. and Marsili, M. (2001). Data clustering and noise undressing of correlation matrices. Phys. Rev. E , 63(6):061101.
- [12] Giada, L. and Marsili, M. (2002). Algorithms of maximum likelihood data clustering with applications. Physica A: Statistical Mechanics and its Applications , 315(3):650 - 664.
- [13] Gopikrishnan, P., Meyer, M., Amaral, L., and Stanley, H. (1998). Inverse cubic law for the distribution of stock price variations. Eur. Phys. J. B , 3(2):139-140.
- [14] Guo, X., Zhang, H., and Tian, T. (2018). Development of stock correlation networks using mutual information and financial big data. PLOS ONE , 13(4):e0195941.
- [15] Hendricks, D., Gebbie, T., and Wilcox, D. (2016a). Detecting intraday financial market states using temporal clustering. Quantitative Finance , 16(11):1657-1678.
- [16] Hendricks, D., Gebbie, T., and Wilcox, D. (2016b). Highspeed detection of emergent market clustering via an unsupervised parallel genetic algorithm. South African Journal of Science , 112(1/2):9.
- [17] Jain, A. K. (2010). Data clustering: 50 years beyond K-means. Pattern Recognition Letters , 31(8):651 - 666.
- [18] Javed, M. A., Younis, M. S., Latif, S., Qadir, J., and Baig, A. (2018). Community detection in networks: A multidisciplinary review. Journal of Network and Computer Applications , 108:87-111.
- [19] Kalayci, C. B., Ertenlice, O., and Akbay, M. A. (2019). A comprehensive review of deterministic models and applications for mean-variance portfolio optimization. Expert Systems with Applications , 125:345-368.
- [20] Khan, K., Rehman, S. U., Aziz, K., Fong, S., and Sarasvady, S. (2014). DBSCAN: Past, present and future. In The Fifth International Conference on the Applications of Digital Information and Web Technologies (ICADIWT 2014) , pages 232-238.
- [21] Kim, J., Lee, W., Song, J. J., and Lee, S.-B. (2017). Optimized combinatorial clustering for stochastic processes. Cluster Comput , 20(2):1135-1148.
- [22] Liao, T. W. (2005). Clustering of time series data-a survey. Pattern Recognition , 38(11):1857-1874.
- [23] Marsili, M. (2002). Dissecting financial markets: sectors and states. Quantitative Finance , 2(4):297-302.
- [24] Marti, G., Andler, S., Nielsen, F., and Donnat, P. (2016). Clustering Financial Time Series: How Long is Enough? In Proceedings of the Twenty-Fifth International Joint Conference on Artificial Intelligence , IJCAI'16, pages 2583-2589. AAAI Press. event-place: New York, New York, USA.
- [25] Marti, G., Nielsen, F., Bi´kowski, M., and Donnat, P. n (2017). A review of two decades of correlations, hierarchies, networks and clustering in financial markets. arXiv preprint arXiv:1703.00485 .
- [26] Mastromatteo, I. and Marsili, M. (2011). On the criticality of inferred models. Journal of Statistical Mechanics: Theory and Experiment , 2011(10):P10012.
- [27] McInnes, L., Healy, J., and Astels, S. (2017). hdbscan: Hierarchical density based clustering. Journal of Open Source Software , 2(11):205.
- [28] McInnes, L., Healy, J., Saul, N., and Groibberger, L. (2018). UMAP: Uniform Manifold Approximation and Projection. The Journal of Open Source Software , 3(29):861.
- [29] Murtagh, F. and Contreras, P. (2012). Algorithms for hierarchical clustering: an overview. WIREs Data Mining and Knowledge Discovery , 2(1):86-97.
- [30] Namaki, A., Shirazi, A. H., Raei, R., and Jafari, G. R. (2011). Network analysis of a financial market based on genuine correlation and threshold method. Physica A: Statistical Mechanics and its Applications , 390(21):38353841.
- [31] Noh, J. D. (2000). Model for correlations in stock markets. Phys. Rev. E , 61(5):5981-5982.
- [32] Ozaki, N., Tezuka, H., and Inaba, M. (2016). A Simple Acceleration Method for the Louvain Algorithm. undefined .
- [33] Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., and Duchesnay, E. (2011). Scikit-learn: Machine Learning in Python. Journal of Machine Learning Research , 12:2825-2830.
- [34] Pollet, J. M. and Wilson, M. (2010). Average correlation and stock market returns. Journal of Financial Economics , 96(3):364-380.
- [35] Prado, M. L. d. (2016). Building Diversified Portfolios that Outperform Out of Sample. The Journal of Portfolio Management , 42(4):59-69.
- [36] Ronhovde, P. and Nussinov, Z. (2010). Local resolutionlimit-free Potts model for community detection. Phys. Rev. E , 81(4):046114.
- [37] Santos, J. M. and Embrechts, M. (2009). On the Use of the Adjusted Rand Index as a Metric for Evaluating Supervised Classification. In ICANN 2009. Lecture Notes in Computer Science, vol 5769. , In: Alippi C., Polycarpou M., Panayiotou C., Ellinas G. (eds) Artificial Neural Networks - ICANN 2009. Springer, Berlin, Heidelberg.
- [38] Steinley, D. (2006). K-means clustering: A half-century synthesis. British Journal of Mathematical and Statistical Psychology , 59(1):1-34.
- [39] Traag, V. A. (2015). Faster unfolding of communities: Speeding up the Louvain algorithm. Phys. Rev. E , 92(3):032801.
- [40] Tumminello, M., Lillo, F., and Mantegna, R. N. (2007). Hierarchically nested factor model from multivariate data. EPL , 78(3):30006.
- [41] Tumminello, M., Lillo, F., and Mantegna, R. N. (2010). Correlation, hierarchies, and networks in financial markets. Journal of Economic Behavior &amp; Organization , 75(1):40-58.
- [42] Wilcox, D. and Gebbie, T. (2007). An analysis of crosscorrelations in an emerging market. Physica A: Statistical Mechanics and its Applications , 375(2):584 - 598.

- [43] Wilcox, D. and Gebbie, T. (2014). Hierarchical Causality in Financial Economics. https://ssrn.com/abstract=2544327 .
- [44] Wu, F. Y. (1982). The Potts model. Rev. Mod. Phys. , 54(1):235-268.
- [45] Yelibi, L. (2020a). Agglomerative Likelihood Clustering. https://github.com/lyelibi/ALC .
- [46] Yelibi, L. (2020b). One Factor Timeseries Model. https://github.com/lyelibi/timeseries gen .
- [47] Yelibi, L. and Gebbie, T. (2020). Fast SuperParamagnetic Clustering. Physica A: Statistical Mechanics and its Applications , 551:124049.
- [48] Zalesky, A., Fornito, A., and Bullmore, E. (2012). On the use of correlation as a measure of network connectivity. NeuroImage , 60(4):2096-2106.
- [49] Zhang, Y., Li, X., and Guo, S. (2018). Portfolio selection problems with Markowitz's mean-variance framework: a review of literature. Fuzzy Optim Decis Making , 17(2):125-158.

## Appendix A: Algorithm: Agglomerative Likelihood Clustering

Here we provide ALC's pseudocode which: 1.) generates clustering candidates, 2.) evaluates the likelihood of the candidate configurations L c [11] using Eqn. (7), and then 3.) selects the best candidates and discards the others.

## TABLE 1

Algorithm 1 Pseudo-code for a ALC implementation (Sec. (III)), and 'alc.py' in [45]

glyph[negationslash]

1: INPUT: Correlation Matrix, OUTPUT: Tracker

2: Produce population of size N

3: while N ≥ 2 do

4: Pick label i in population

5: for j = i in population

do

6: create label k by merging labels i and j

7: Pick k by optimizing for ∆ L c

8: add k, remove i and j from population.

9: end for

10:

Stop if ∆ L c ≤ 0 and output Tracker

11:

for

j in population do

12:

store correlation C kj in correlation matrix

13:

update cluster composition in Tracker

14: end for

15: end while

## Appendix B: Data Generative Process: Synthetic Correlated Time-Series

The Noh Ansatz [11, 31] offers a powerful stochastic process model for clustering purposes. We make use of the following equation (equivalent to Eqn. (2) ) as a way of generating correlated time-series:

$$\bar { \xi } _ { i } ( d ) = \frac { \sqrt { g _ { s _ { i } } } \eta _ { s _ { i } } ( d ) + \epsilon _ { i } ( d ) } { \sqrt { 1 + g _ { s _ { i } } } } \, \quad \quad ( B 1 )$$

The process is described in Table (2). In Fig. (2a) We plotted 3000 simulated time-series of D ∈ (20 60 250) , , with 10 clusters each of size 300. Time-series visualizations lack interpretability especially given apparent noise and chaos. We use additional tools such as Minimum Spanning Trees (MST) and dimensionality reduction methods such as UMAP to provide better visualization. Shown in Fig. (2b) and (2c) are respectively the MSTs for the true correlation matrix, and the estimated of the time-series returns in Fig. (2a). Because MSTs rely on nearest neighbors to build the graph, here the input is d ij = 1 -ρ ij . The MSTs are clearly capable of dissociating clusters, and one could, in this simple case, discern all 10 clusters without much prior information. A minor remark is in the topology of the MST which differs from the true to the estimated correlation matrix: estimating correlations from data introduces noise which is reflected in its topology 12 .

- (2)1 Cluster Number : Define values for number of cluster C , and size of clusters s and obtain N = s ∗ C the number of time-series in the data-set. Pick a time-series length D
- (2)2 Spin-Labels : Create a list of array of spin-labels with the C labels
- (2)3 Random Effects : Create a CxD array η ∼ N (0 , 1), and another NxD array glyph[epsilon1] ∼ N (0 , 1). η and glyph[epsilon1] respectively capture the daily cluster and individual object random effects.
- (2)4 Fix Intra-cluster Binding Strength : Pick a value for g s per cluster: it is not needed that g s be identical for every clusters. As a matter of fact, real noisy systems will display various g s values. However fixing g s simplifies the process and increases interpretability.
- (2)5 Compute Returns : Create a NxD array ξ and compute the daily returns using Eqn. (B1) by looping over the clusters, and the time-series within the clusters.

TABLE 2: Implementation of Noh Ansatz model of correlated time-series (See code at [46] )

Using UMAP we plot the same data in Fig. (2a) but projected on a 2D manifold in Fig. (8) for visual purposes

12 Noise due to the estimation of the correlation matrix are typically cleaned to a certain extent using Random Matrix Theory methods [5, 42]

FIG. 8: From top to bottom g s is varied from 0.1 to 0.3 (See model in Appendix (B) ) and from left to right, time-series length from 20 to 250. The result is the UMAP [28] 2D embed of 3000 correlated time-series sorted within 10 clusters.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00951v4-figure-c894b50f-3c4e-40c6-9622-950e5362912f.jpg)

only. We expect UMAP to capture the local structure of the data: From left to right the time-series lengths are 20 to 250, while from top to bottom g s values range from 0.1 to 0.3. The noise level is captured by g s : clusters with g s → 0 are diffuse and close to each other whereas those with g s → 1 show increasingly high density. Similarly short time-series induce spurious correlations and noise within clusters. The data-sets remain modular and very diffuse but as time-series length increases they spread out and thus should become easier to separate by clustering algorithms.
Lionel Yelibi Tim Gebbie
Agglomerative Likelihood Clustering
We consider the problem of fast time-series data clustering. Building on
previous work modeling the correlation-based Hamiltonian of spin variables we
present an updated fast non-expensive Agglomerative Likelihood Clustering
algorithm (ALC). The method replaces the optimized genetic algorithm based
approach (f-SPC) with an agglomerative recursive merging framework inspired by
previous work in Econophysics and Community Detection. The method is tested on
noisy synthetic correlated time-series data-sets with built-in cluster
structure to demonstrate that the algorithm produces meaningful non-trivial
results. We apply it to time-series data-sets as large as 20,000 assets and we
argue that ALC can reduce compute time costs and resource usage cost for large
scale clustering for time-series applications while being serialized, and hence
has no obvious parallelization requirement. The algorithm can be an effective
choice for state-detection for online learning in a fast non-linear data
environment because the algorithm requires no prior information about the
number of clusters.