Graphical Models in Heavy-Tailed Markets
Year: 2021

Abstract:
Heavy-tailed statistical distributions have long been considered a more realistic statistical model for the data generating process in financial markets in comparison to their Gaussian counterpart. Nonetheless, mathematical nuisances, including nonconvexities, involved in estimating graphs in heavy-tailed settings pose a significant challenge to the practical design of algorithms for graph learning. In this work, we present graph learning estimators based on the Markov random field framework that assume a Student-t data generating process. We design scalable numerical algorithms, via the alternating direction method of multipliers, to learn both connected and k-component graphs along with their theoretical convergence guarantees. The proposed methods outperform state-of-the-art benchmarks in an extensive series of practical experiments with publicly available data from the S&P500 index, foreign exchanges, and cryptocurrencies.

heading: 1 Introduction
text: Graph learning frameworks are often designed based on the assumption that the observed graph signals are Gaussian distributed [1–9]. While such assumption for graphical models has found great success in many practical areas, which includes brain network analysis [10], psychological networks [11], and single-cell sequencing [12], it inherently neglects scenarios where there may exist outliers or the underlying data is naturally heavy-tailed distributed. As a consequence, those methods often lack robustness and may not succeed in capturing a meaningful representation of the underlying graph [13].
Data from financial instruments are well-known examples of such scenarios where heavy-tailedness and skewedness are present [14–19]. In addition, there has been a growing interest in methods for estimating graphical models in financial markets, which hence demands the development of scalable and robust learning algorithms [20].
Perhaps one of the most prominent applications, clustering financial time-series via graph techniques has been an active research topic [20–24]. Nonetheless, current techniques rely on the assumption that the underlying graph has a tree structure, which does bring advantages due to its hierarchical clustering properties, but also have been shown to be unstable [25–27] and not suitable when the data is not Gaussian distributed [28].
Motivated by practical challenging applications in finance, such as clustering of financial instruments and network estimation, we investigate the problem of learning graph matrices whose structure follow that of a Laplacian matrix of an undirected weighted graph for which the data generating process is assumed to be Student-t distributed. In particular, the main contributions of this paper are as follows:
• We propose a novel formulation for learning undirected weighted graphs under the assumption that the data generating process is Student-t distributed. We solve the underlying
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
learning problem via a carefully designed numerical algorithm based on the alternating direction method of multipliers (ADMM), along with the establishment of its theoretical convergence guarantees. We note that the proposed algorithm can be easily extended to account for additional linear constraints on the graph weights.
• We extend the proposed framework to account for heavy-tails and k-component graphs simultaneously, which enables a novel method for clustering financial time-series.
• We present extensive practical results, with real-world data from the US stock market, foreign exchanges, and cryptocurrencies, that showcase clear advantages of including heavy-tail assumptions into graph learning frameworks when compared to state-of-the-art, Gaussian-based methods.
Notation: Matrices (vectors) are denoted by bold, italic, capital (lowercase) roman letters likeX , x. Vectors are assumed to be column vectors. The (i, j) element of a matrixX ∈ Rn×p is denoted as Xij . The i-th element of a vector x is denoted as xi. The i-th row of X is denoted as xi ∈ Rp×1. Given a symmetric matrixA, λi(A) and λmax(A) denote the i-th smallest and maximum eigenvalue of A, respectively. The Moore-Penrose inverse of A is denoted as A†. The Frobenius norm of a matrix A is denoted as ‖A‖F = √ tr (A>A). The operator Diag : Rp → Rp×p creates a diagonal matrix with the elements of an input vector along its diagonal. The operator diag : Rp×p → Rp extracts the diagonal of a square matrix. For x ∈ Rp, ‖x‖∞ = maxi|xi|. (x)+ denotes the projection on to the nonnegative orthant, i.e., (x)+ = max (0,x).

heading: 2 Background & Related Works
text: An undirected, weighted graph is denoted as a triple G = (V, E ,W ), where V = {1, 2, . . . , p} is the node set, E ⊆ {{u, v} : u, v ∈ V, u 6= v} is the edge set, that is, a subset of the set of all possible unordered pairs of nodes such that {u, v} ∈ E iff there exists a link between nodes u and v. W ∈ Rp×p+ is the symmetric weighted adjacency matrix that satisfies Wii = 0,Wij > 0 iff {i, j} ∈ E and Wij = 0, otherwise. The combinatorial, unnormalized graph Laplacian matrix L is defined, as L ,D −W , whereD , Diag(W1) is the degree matrix.
A p-dimensional, real-valued, Gaussian random variable x, with mean vector E [x] , µ and rankdeficient precision matrix L, is said to form a Laplacian constrained Gaussian Markov random field (LGMRF) [9, 29–31] of rank p− k, k ≥ 1, with respect to a graph G, when its probability density function is given as
p(x) ∝ √ det∗ (L) exp { −1
2 (x− µ)>L(x− µ)
} , (1)
where det∗(L) is the pseudo-determinant of L, i.e., the product of its positive eigenvalues [32].
Assume we are given n observations of x, i.e., X = [x1,x2, . . . ,xn] >, X ∈ Rn×p, xi ∈ Rp×1. The goal of graph learning algorithms is to learn a Laplacian matrix, or equivalently an adjacency matrix, given only the data matrixX , i.e., often without any knowledge of E . To that end, the penalized maximum likelihood estimator (MLE) of the Laplacian constrained precision matrix of x, on the basis of the observed dataX , is:
minimize L 0 tr (LS)− log det∗ (L) + h(L), subject to L1 = 0, Lij = Lji ≤ 0,
(2)
where S is a similarity matrix, e.g., the sample covariance (or correlation) matrix S ∝X>X , and h is a regularization function to promote certain properties on L such as sparsity or low-rankness.
Even though Problem (2) is convex, provided we assume a convex choice for h, it is not adequate to be solved by disciplined convex programming languages, such as cvxpy [33], particularly due to scalability issues related to the computation of the term log det∗(L) [6, 7]. Indeed, recently, considerable efforts have been directed towards the design of scalable, iterative algorithms based on block coordinate descent [34], majorization-minimization (MM) [35, 36], and ADMM [37] to solve Problem (2) in an efficient fashion, e.g., [6] and [7].
Estimators based on Gaussian assumptions have been proposed for connected graphs [2, 4–7]. Some of their properties, such as sparsity, are yet being investigated [9, 31]. The authors in [29] and [38] proposed optimization programs for learning the class of k-component graphs, as such class is an appealing model for clustering tasks due to the spectral properties of the Laplacian matrix. However, a major shortcoming in their formulations is the lack of constraints on the degrees of the nodes, which allows for trivial solutions, i.e., graphs with isolated nodes.
Elliptical losses along with linear structural constraints that retain the positive-definiteness of the estimated covariance matrix have been proposed in the literature [39]. In this work, however, we address the case of Laplacian constraints, which lead to positive-semidefinite precision matrices, and nonconvex k-component structural constraints.

heading: 3 Proposed Formulations & Algorithms
text: In this section, we propose optimization formulations and an iterative algorithm to learn a Laplacian matrix from heavy-tailed assumptions. With that goal, we express the Laplacian matrix via its linear operator, i.e., L = Lw [29], wherew ∈ Rp(p−1)/2 is the vectorized form of the upper triangular part of the adjacency matrix, also known as the vector of graph weights. In addition, we use the fact that, for connected graphs, it follows that det∗(Lw) = det(Lw + J), J , 1p11 > [6].
In order to address the inherent heavy-tailed nature of financial market data [40], we consider the Student-t distribution under the improper Markov random field assumption [41] with Laplacian structural constraints, that is, we assume the data generating process to be modeled as multivariate zero-mean Student-t distribution, whose probability density function can be written as
p(x) ∝ √ det∗(Θ) ( 1 + x>Θx
ν
)−ν + p 2
, ν > 2, (3)
where Θ is a positive-semidefinite inverse scatter matrix modeled as a combinatorial graph Laplacian matrix and ν is the number of degrees of freedom, which measures the rate of decay of the tails.
This results in a robustified version of the MLE for connected graph learning, i.e.,
minimize w≥0,Θ 0
p+ ν
n n∑ i=1 log ( 1 + x>i Lwxi ν ) − log det (Θ + J) ,
subject to Θ = Lw, dw = d, (4)
where d : Rp(p−1)/2 → Rp is the degree operator defined as dw , diag (Lw). The constraint dw = d enables the learning of additional graph structures such as regular graphs and it is crucial for k-component graphs, as discussed in Section 3.2.
From a theoretical perspective, the Student-t model naturally yields sparse graphs. Comparing the objective function in Problem (4) to that of Problem (2), we note that the Student-t contains a log(·) term in place of a linear term of the graph weights. The usage of a log function to promote sparsity is closely related to the iteratively reweighted `1-norm as an approximation for the `0-norm problem [42]. Problem (4) is, in general, nonconvex due to the summation of log terms and hence it is challenging to be considered directly. Hence, we design an iterative algorithm based on the ADMM framework.

heading: 3.1 ADMM Solution
text: The partial augmented Lagrangian function of Problem (4) is given as
Lρ(Θ,w,Y ,y) = p+ ν
n n∑ i=1 log ( 1 + x>i Lwxi ν ) − log det (Θ + J) + 〈y, dw − d〉+ ρ 2 ‖dw − d‖22 + 〈Y ,Θ− Lw〉+ ρ 2 ‖Θ− Lw‖2F , (5)
where Y and y are the dual variables associated with the equality constraints Θ = Lw and dw = d, respectively. Note that we deal with the constraints w ≥ 0 and Θ 0 directly, hence there are no dual variables associated with them.
Given wl and Y l, the subproblem for Θ can be written as
Θl+1 = arg min Θ 0 − log det (Θ + J) + 〈Θ,Y l〉+ ρ 2
∥∥Θ− Lwl∥∥2 F , (6)
whose closed-form solution is given by Lemma 1.
Lemma 1 The global minimizer of problem (6) is [43, 44]
Θl+1 = 1 2ρ U ( Γ + √ Γ2 + 4ρI ) U> − J , (7)
where UΓU> is the eigenvalue decomposition of ρ ( Lwl + J ) − Y l.
Given Θl+1, Y l, and yl, the subproblem for w can be formulated as
minimize w≥0
ρ 2 w> (d∗d + L∗L)w −
〈 w,L∗ ( Y l + ρΘl+1 ) − d∗ ( yl − ρd )〉 + p+ ν
n n∑ i=1 log ( 1 + x>i Lwxi ν ) , (8)
where d∗ and L∗ are the adjoint operators of the degree and Laplacian operators, respectively. In general, subproblem (8) is nonconvex due to the concave nature of the logarithm function. Hence, we resort to the MM method [36] to find a stationary point of subproblem (8). We proceed by constructing a global upper-bound of the objective function of (8) at point wj ∈ Rp(p−1)/2+ as
g(w,wj) = g(wj ,wj) + 〈w −wj ,∇wf(wj)〉+ µ
2
∥∥w −wj∥∥2 2 , (9)
where f is the objective function in the minimization in (8), its gradient is given as ∇wf(wj) = aj + bj , where
aj = L∗ ( S̃j − Y l − ρ ( Θl+1 − Lwj )) , (10)
bj = d∗ ( yl − ρ ( d− dwj )) , (11)
where S̃j , 1n ∑n i=1
(p+ ν)xix > i
〈wj ,L∗(xix>i )〉+ ν is a weighted sample covariance matrix, and µ =
ρλmax (d ∗d + L∗L), and the maximum eigenvalue of d∗d + L∗L is given by Lemma 2, whose proof is presented in the Supplementary Material.
Lemma 2 The maximum eigenvalue of the matrix d∗d + L∗L is given as λmax (d ∗d + L∗L) = 2(2p− 1). (12)
The vector of graph weights w can then be updated by minimizing the function g constructed in (9), which is tantamount to solving the following nonnegative, quadratic-constrained, strictly convex problem:
wj+1 = arg min w≥0
ρ(2p− 1) ∥∥w −wj∥∥2
2 + 〈w,aj + bj〉, (13)
whose unique solution can be readily obtained via its KKT optimality conditions and is given as
wj+1 = ( wj − a j + bj
2ρ(2p− 1)
)+ , (14)
that is, a projected gradient descent step with learning rate (2ρ(2p− 1))−1. Thus, we iterate (14) in order to obtain a stationary point, wl+1, of Problem (8). In practice, we observe that a few iterations are sufficient to retrieve wl+1.
The dual variables Y and y are updated via gradient ascent steps. Algorithm 1 summarizes the implementation to find a stationary point of Problem (4). We present Theorem 3, proved in the Supplementary Material, which establishes the convergence of Algorithm 1. Theorem 3 Algorithm 1 converges subsequently for any sufficiently large ρ, that is, the sequence{( Θl,wl,Y l,yl )} generated by Algorithm 1 has at least one limit point, and each limit point is a stationary point of Problem (4).
Algorithm 1: Student-t Graph Learning Data: Data matrixX ∈ Rn×p, initial estimate of the graph weights w0, desired degree vector d, penalty parameter ρ > 0, degrees of freedom ν > 2, convergence tolerance > 0 Result: Graph Laplacian estimation: Lw?
1 initialize Y = 0, y = 0 2 l← 0 3 while ∥∥rl∥∥∞ > or ∥∥sl∥∥∞ > do 4 . update Θl+1 via (7) 5 . update wl+1 by iterating (14) 6 . update Y l+1 = Y l + ρ ( Θl+1 − Lwl+1
) 7 . update yl+1 = yl + ρ ( dwl+1 − d
) 8 . compute residual rl+1 = Θl+1 − Lwl+1 9 . compute residual sl+1 = dwl+1 − d
10 l← l + 1 11 end
3.2 An extension to k-component graphs
The graph learning formulation proposed in (4) is applicable to learn connected graphs. Learning graphs with k components, k assumed to be known, poses a considerably higher challenge, as the dimension of the nullspace of the Laplacian matrix Lw is equal to the number of components of the graph [45]. One way to achieve this requirement is by imposing the constraint rank (Lw) = p− k, which is nonconvex and nondifferentiable, in the maximum likelihood problem generated by (3). We instead relax this rank constraint by noting that via Fan’s theorem [46], we have
k∑ i=1 λi (Lw) = minimize V ∈Rp×k,V >V =I tr ( V >LwV ) . (15)
Thus, by using the right hand side of (15) as a regularization term, we are able formulate the following optimization problem to learn a Student-t k-component graph:
minimize w≥0,Θ 0,V
p+ ν
n n∑ i=1 log ( 1 + x>i Lwxi ν ) − log det∗ (Θ) + ηtr(LwV V >),
subject to Θ = Lw, rank(Θ) = p− k, dw = d, V >V = I, V ∈ Rp×k. (16)
The partial augmented Lagrangian function of Problem (16) can be expressed as
Lρ(Θ,w,V ,Y ,y) = p+ ν
n n∑ i=1 log ( 1 + x>i Lwxi ν ) − log det∗ (Θ) + ηtr ( LwV V > ) + 〈y, dw − d〉+ ρ
2 ‖dw − d‖22 + 〈Y ,Θ− Lw〉+
ρ 2 ‖Θ− Lw‖2F . (17)
Given wl and Y l, the subproblem for Θ can be written as
Θl+1 = arg min rank(Θ)=p−k − log det∗(Θ) + 〈Θ,Y l〉+ ρ 2
∥∥Θ− Lwl∥∥2 F , (18)
which is nearly the same as (6). Its solution is obtained as
Θl+1 = 1 2ρ U ( Γ + √ Γ2 + 4ρI ) U>, (19)
except that now UΓU> is the eigenvalue decomposition of ρLwl − Y l, with Γ having the largest p− k eigenvalues along its diagonal and U ∈ Rp×(p−k) contains the corresponding eigenvectors. The subproblem to obtainwl+1 is virtually the same as in (8) except for the additional linear term ηtr(LwV lV l>). Hence, its update is also a projected gradient descent step, alike (14) where
aj , L∗ ( S̃j + ηV lV l> − Y l − ρ ( Θl+1 − Lwj )) . (20)
Given wl+1, we have the following subproblem for V :
minimize V ∈Rp×k,V >V =I
tr ( V >Lwl+1V ) , (21)
whose closed-form solution is given by the k eigenvectors associated with the k smallest eigenvalues of Lwl+1 [47, 48]. Algorithm 2 summarizes the implementation to find a stationary point of Problem (16), and its convergence is established through Theorem 4, whose proof is presented in the Supplementary Material. Theorem 4 Algorithm 2 converges subsequently for any sufficiently large ρ, that is, the sequence{( Θl,wl,V l,Y l,yl )} generated by Algorithm 2 has at least one limit point, and each limit point is a stationary point of Problem (16).
Algorithm 2: k-component Student-t graph learning Data: Data matrixX ∈ Rn×p, initial estimate of the graph weights w0, number of graph
components k, desired degree vector d, degrees of freedom ν, rank hyperparameter η > 0, penalty parameter ρ > 0, tolerance > 0
Result: Laplacian estimation: Lw? 1 initialize Y = 0, y = 0 2 l← 0 3 while
∥∥rl∥∥∞ > or ∥∥sl∥∥∞ > do 4 . update Θl+1 via (19) 5 . update wl+1 as in (14) with aj given in (20) 6 . update V l+1 as in (21) 7 . update Y l+1 = Y l + ρ ( Θl+1 − Lwl+1
) 8 . update yl+1 = yl + ρ ( dwl+1 − d
) 9 . compute residual rl+1 = Θl+1 − Lwl+1
10 . compute residual sl+1 = dwl+1 − d 11 l← l + 1 12 end

heading: 4 Experiments
text: To evaluate the performance of the proposed graph learning algorithms, we perform experiments using historical daily price time series data, available in Yahoo! FinanceTM, from financial instruments in three scenarios: (i) stocks belonging to the S&P500 index, (ii) foreign exchange markets, and (iii) cryptocurrencies. We start by constructing the log-returns data matrix, i.e., a matrixX ∈ Rn×p, where n is the number of log-return observations and p is the number of instruments, as
Xi,j = logPi,j − logPi−1,j , (22)
where Pi,j is the closing price of the j-th instrument at the i-th day.
Benchmarks: We compare our proposed algorithms with state-of-the-art, Gaussian distribution-based methods for connected graphs, namely GLE [7] and NGL [9], which use `1-norm and minimax concave penalty regularizations, respectively; and CLR [38] and SGL [29] that consider k-component graphs. For a fair comparison among algorithms, we set the degree vector d equal to 1 for the proposed algorithms, i.e., we do not consider any prior information on the degree of nodes. In our ADMM algorithms, we set the penalty parameter to ρ = 1 and the hyperparameter η in (16) is adaptively increased until the rank constraint is satisfied. For GLE and NGL, we use grid search on the sparsity hyperparameter such that the resulting graph yields the highest modularity value. The graph weights in Algorithm 1 and 2 are initialized using the same procedure as in [8].
Our goal with the experiments that follow is to verify whether the heavy-tail assumption provides an improved version of the learned graph, which is evaluated based on the modularity1 of the estimated
1Modularity measures the strength of separability of a graph into groups [49].
graph and the graph visualization. In addition, for the task of clustering stocks, we analyze whether the learned graphs agree with industry standards of sector classification set by the Global Industry Classification Standard (GICS) [50, 51]. 2

heading: 4.1 Communities in S&P500 Stocks
text: In this experiment, we consider S&P500 stocks belonging to three sectors, namely, Communication Services (red), Utilities (blue), and Real Estate (green), totalling p = 82 stocks, during the time horizon from Jan. 3rd 2014 to Dec. 29th 2017, resulting in n = 1006 observations. In order to obtain descriptive insights on this dataset, we measure its degree of heavy-tailedness and annualized volatility3. The former is obtained by fitting the degrees of freedom of a Student-t distribution to the matrix of log-returns, whereby we obtain ν ≈ 5.5 and σ ≈ 21%. This scenario can be considered as having a moderate amount of heavy-tailedness.
Figure 1 depicts the learned connected graphs on the aforementioned time periods. It can be readily noticed that the graph learned with the Student-t distribution (Figure 1c) is sparser than those learned with the Gaussian assumption (Figure 1a and 1b), which results from the fact that the Gaussian distribution is more sensitive to outliers. Moreover, the Student-t graph presents a higher degree of interpretability as measured by its modularity value. In addition, a larger number of inter-sector connections, as indicated by gray-colored edges, which are often spurious from a practical perspective, are present in the graphs learned by NGL and GLE. Sparsity regularization provides a means to remove edges between nodes in the presence of data with outliers and possibly increasing the modularity of the resulting graph. However, they bring the additional task of tunning hyperparameters, which is often repetitive and impractical for real-time applications. A cleaner graph, without the need for postprocessing or additional regularization, is obtained directly by using the Student-t assumption.

heading: CNP
text: Figure 2 illustrates the learned 3-component graphs during the time from Jan. 3rd 2014 until Dec. 29th 2017. We can notice that SGL (Figure 2a) and CLR (Figure 2b) are unable to separate the stocks in a way that agrees with their sector information as given by GICS. In addition, the high number of spurious connections (gray-colored edges) are uncharacteristic of the actual expected behavior in stock markets. Figure 2c displays the graph learned by the proposed Algorithm 2, where it presents not only a higher modularity value, but also a sparser, more plausible representation of an actual network of stocks with three sectors.

heading: 4.1.1 Impact of COVID-19 in the US Stock Market
text: We collect price data of p = 85 stocks belonging to three sectors, namely, Communication Services (red), Utilities (blue), and Real Estate (green) from Apr. 22nd 2019 to Dec. 30th 2020, resulting in n = 429 observations. The degrees of freedom and annualized volatility during this period, which includes the financial crisis caused by the COVID-19 pandemic, were ν ≈ 2.89 and σ ≈ 41%, indicating a strongly heavy-tailed scenario.
2More often than not, stocks have impacts on multiple sectors, e.g., the evident case of technology companies, whose influence affect the prices of stocks not only in their own sector, but spans across multiple sectors. Therefore, GICS or other sector classification systems cannot be considered as absolute ground-truth labels.
3The annualized volatility is computed as σ = √
252 p ∑p i=1 σi, where σi is the daily sample standard
deviation of the i-th stock.
In Figure 3, we observe that the modularity value of the Gaussian-based graphs are severely reduced as compared to the results presented in Figure 1, while the opposite occurs with the Student-t graph. That may be explained by the increase in variability and outliers in the data during the financial crisis caused by the COVID-19 pandemic, which is better modelled by a Student-t distribution. In addition, the learned graph in Figure 3c clearly displays expected behaviors such as the strong correlation between pairs of companies including {Zoom Video Communications Inc., Netflix Inc.}, and {Twitter Inc., Facebook Inc.}, which are not as evident in Figures 3a and 3b.

heading: EQIX

heading: SBAC

heading: 4.2 Communities in Foreign Exchange Markets
text: We query foreign exchange data from the 34 most traded currencies between the period from Jan. 2nd 2019 to Dec. 31st 2020, totalling n = 522 observations. The data matrix is composed by the log-returns of the currencies prices with respect to the United States Dollar. Similar to the previous experiment, we compute the degrees of freedom of a Student-t distribution fitted to the log-returns data matrix, whereby we obtain ν ≈ 4.6, which represents a scenario with considerable amount of heavy-tailedness. Unlike in the experiment involving S&P500 stocks, there are no classification standard for currencies, hence we rely on a community detection algorithm [52] in order to create classes within the learned graph. In particular, the algorithm in [52] takes as input the learned Laplacian matrix of the graph and outputs a membership assignment that maximizes the modularity of the graph.
Figure 4 displays the learned graphs. As it can be observed, the Student-t graph (Figure 4c) is sparser, more interpretable, and has a higher modularity value than that of the Gaussian-based graphs (Figure 4a and 4b). In addition, the expected correlation between currencies of locations geographically close to each other, e.g., {Hong Kong SAR, China}, {Taiwan, South Korea}, and {Poland, Czech Republic} are significantly more evident for the Student-t graph. Figure 5 depicts the learned 9-component graphs of currencies during the time window from Jan. 2nd 2019 to Dec. 31st 2020. It can be observed that the graph learned by the proposed Algorithm 2 (Figure 5c) presents a finer structure and a higher modularity value than those learned by SGL (Figure 5a) and CLR (Figure 5b). In addition, the learned graph in Figure 5c presents more reasonable clusters such as {New Zealand, Australia} and {Poland, Czech Republic, Hungary}, which are not
separated in the Gaussian-based graphs. More critically, we can observe that SGL and CLR allow the existence of isolated nodes in the learned graphs. In our proposed algorithm, we avoid such solutions by imposing linear constraints on the degree of the graph.

heading: 4.3 Communities in Cryptocurrencies
text: We query daily prices of the p = 41 most traded cryptocurrencies during the period starting from Aug. 1st 2017 to Dec. 1st 2020, which amounts to n = 1218 observations. The degrees of freedom during this time frame was measured as ν ≈ 3, which is tantamount to a strong heavy-tail scenario. Figure 6 shows the learned graphs by GLE, NGL, and our proposed Algorithm 1. While the graphs in Figures 6a and 6b present a small modularity value and contain a large number of edges, which impairs interpretability, the resulting proposed graph in Figure 6c reveals a refined representation of the interactions between pairs of cryptocurrencies, which is possibly more aligned with the actual market scenario. As an example, the link between Bitcoin (BTC) and Litecoin (LTC), a Bitcoin spinoff established in 2011, is substantially more evident in Figure 6c.

heading: XEM
text: Figure 7 shows the learned 7-component graphs of cryptocurrencies during the aforementioned time window. As in the previous experiments with foreign exchange data, SGL shows isolated nodes in the learned graph (Figure 7a) and CLR contains a large number of spurious connections in the main cluster (Figure 7b), whereas the graph learned via Algorithm 2 (Figure 7c) has the largest modularity value. Interestingly, while all three methods agree to cluster {Dogecoin (DOGE), Verge (XVG), Siacoin (SC), DigiByteCoin (DGB)}, which may be related to their similar initial release dates, only our proposed algorithm clusters together the coins that mainly focus on privacy and anonymity features, i.e., {Monero (XMR), Zcash (ZEC), DASH}.

heading: BNT

heading: MONA

heading: 4.4 Effect of Algorithm Initialization
text: Initialization is a critical phase in any iterative algorithm, especially when dealing with nonconvex optimization problems involving nonconvex constraints as it is the case of Algorithm 2. To verify the stability of Algorithm 2, we initialize it as a fully connected graph with equal weights, i.e., w0 ∝ 1. Figure 8 illustrates that there is no significant difference from the estimated graphs with uniform initial point and the ones reported in the previous section. This result provides evidence that Algorithm 2 is robust against inaccurate initializations.

heading: 5 Conclusions
text: Heavy-tails are prevalent in time-series of financial markets. Yet, they have been little explored in the context of Laplacian graphical models. In this paper, we have proposed optimization programs to learn graphical models with Laplacian constraints assuming that the data generating process is Student-t distributed. The formulations follow a maximum likelihood approach of a Markov random field, for which we designed ADMM algorithms that converge to a stationary point of the resulting nonconvex problems. The proposed algorithms showed significant gains, measured via the modularity values of the estimated graphs, when compared to state-of-the-art counterparts in real-world scenarios that involved data from the US stock market, foreign exchange markets, and cryptocurrencies.

heading: Acknowledgments
text: We would like to thank the anonymous reviewers for their helpful comments. This work was supported by the Hong Kong GRF 16207820 research grant.

Review summaries:
The paper proposes a new estimator for learning graphical models under the t-distribution assumption for data generation process. For the new estimator, an ADMM-based algorithm is designed to solve the associated optimization problem in order to find the graph structures. The new estimator is evaluated extensively on real-world datasets to demonstrate its effectiveness. The paper also establishes the convergence of the proposed algorithm. 

The paper proposes a framework for learning a graphical model when the data generating distribution is heavy-tailed. More specifically, the author(s) propose a graph learning estimator under Markov Random Field framework when the data is generated from a Student's t-distribution. A scalable algorithm is proposed to handle connected and K-component graphs. The methodology is particularly applicable to financial time series data. The algorithm is shown to be superior in performance compared to existing methods in both simulated and real data applications.

- This paper proposes to use the student-t distribution for data generation in financial markets, which can capture the effect of outliers better than the commonly used Gaussian distribution.   - In particular, the data is modeled using a Laplacian Markov random field with student-t distribution, and a practical ADMM based algorithm to learn the underlying graph structure is proposed.  - The algorithms are designed under two settings namely, connected graphs and k-component graphs, and convergence is proven.   - The proposed approach is shown to outperform the state-of-the-art Gaussian assumption-based approaches with real-world data from stock markets, foreign exchanges, and cryptocurrencies.