Multiscale Causal Structure Learning
Authors: Gabriele D'Acunto
Paolo Di Lorenzo
Sergio Barbarossa
Venue: Accepted by TMLR ()

Abstract: Causal structure learning methods are vital for unveiling causal relationships embedded into observed data. However, the state of the art suffers a major limitation: it assumes that causal interactions occur only at the frequency at which data is observed. To address this limitation, this paper proposes a method that allows structural learning of linear causal relationships occurring at different time scales. Specifically, we explicitly take into account instantaneous and lagged inter-relations between multiple time series, represented at different scales, hinging on wavelet transform. We cast the problem as the learning of a multiscale causal graph having sparse structure and dagness constraints, enforcing causality through directed and acyclic topology. To solve the resulting (non-convex) formulation, we propose an algorithm termed MS-CASTLE, which exhibits consistent performance across different noise distributions and wavelet choices. We also propose a single-scale version of our algorithm, SS-CASTLE, which outperforms existing methods in computational efficiency, performance, and robustness on synthetic data. Finally, we apply the proposed approach to learn the multiscale causal structure of the risk of 15 global equity markets, during covid-19 pandemic, illustrating the importance of multiscale analysis to reveal useful interactions at different time resolutions. Financial investors can leverage our approach to manage risk within equity portfolios from a causal perspective, tailored to their investment horizon.

Published in Transactions on Machine Learning Research (10/2023)
Multiscale Causal Structure Learning
Gabriele D' D'Acunto DIAG, Sapienza University of Rome Centai Institute, Turin, Italy
gabricde.duccnloo@unnnnnn...
Paolo Di Lorenzo
polb..cllovesso@@nnnrrooa...
DIET, Sapienza University of Rome
Sergio Barbarossa
sergio. Stanobrrrrrasso@@iroooo....
DIET, Sapienza University of Rome
Reviewed on OpenReview: https:
// openreview. net/ forum? id= Ub6XILEF9a
Abstract
Causal structure learning methods are vital for unveiling causal relationships embedded into observed data. However, the state of the art suffers a major limitation: it assumes that causal interactions occur only at the frequency at which data is observed. To address this limitation, this paper proposes a method that allows structural learning of linear causal relationships occurring at different time scales. Specifically, we explicitly take into account instantancous and lagged inter-rclations between multiple time series, represented differ- ent scales, hinging on wavelet transform. We cast the problem as the learning of a multiscale causal graph having sporse structure and dogness constraints, enforeing causality through directed and acyelic topology. To solve the resulting (non-convex) formulation, propose an algorithm termed MS-CASTLE, which exhibits consistent performance across differ- ent. noise distributions and wavelet choices. We also propose a single-scale version of our algorithm, SS-CASTLE, which outperforms existing methods in computational efficiency, performance, and robustness on syuthetic data. Finally, we apply the proposed approach to learn the multiscale causal structure of the risk 15 global equity markets, during covid- 19 pandemic, illustrating the importance of multiscale analysis to reveal useful interactions at different time resolutions. Financial investors can leverage our approach to manage risk within equity portfolios from a causal perspective, tailored to their investment horizon.
1 Introduction
The study of causal relationships plays fundamental role in our understanding of complex systems. However, learning causal relationships is a challenging task, and often not even possible to directly act on the systems of interest (e.g., social networks) because of ethics, feasibility, or cost issues. Thus, the ability unravel causal structures from the observed data, also known as causal structure learning, an attractive technology that has received growing attention in the last years, also thanks to the ever-increasing volume of available data, see e.g., [Pearl 2009 Peters et al. 2017 Glymour et al. 2019 Schölkopf et al. 202 In the literature, several works hinged on directed acyclic graphs (DAGs) to represent causal dependencies (i.e., directed edges) between the constituents (e.g., nodes) of the considered system (Vowels el al., 2022). Indeed, the acyclicity requirement represents necessary condition in order to set causes apart from effects; result that cannot be accomplished the presence of feodback loops among the nodes. In case the nodes of the DAG are associated with time series observations, the dependencies represented by the DAG refer causal interactions occurring between the values of the time series within the same timestamp and across

---

Published in Transactions on Machine Learning Research (10/2023)
different time stamps. The formers are called instantancous, whereas the latters, since we can only observe causal interactions coming from the past, are named lagged.
Related works. Causal structure learning algorithms can be classified in accordance with the approach used to infer the associated DAG. In particular, we can identify three main different classes: (i) constraint- based approaches, which run conditional independence tests to validate the presence of an edge between two variables (Spirtes et al.) 2000; Huang et al., 2020g; (ii) score-based methods, which measure the goodness of fit. of graphs according to a given criterion and then use search procedures to explore the solution space Chickering [2002] Huang et al., [0188;; (iii) functional methods, which model a variable in terms of function of its parents (Shimizu et al. 2006 Hoyer et al.] 2008 Hyvärinen et al.,
(Heckerman et al., 1995
2010
Peters et al.) 2014 Bühlmann et al. 2014). Furthermore, recent important contribution came by
reformulating the problem of learning DAG using a suitable continuous non-convex penalty (Zheng et al.
2018
2020).
This enabled the usage of gradient-based methods (Yu et al. 2019; Lachapelle et al. 2020 Ng
ct al. 2020 and reinforcement learning (Zhu et al. 2020) in causal discovery problems.
Whenever we deal with causal inference for time series analysis, we need to take into account time ordering as well. Considering linear models, this leads to the formulation of the structural vector autoregressive model (SVARM) (Kilian & Lütkepohl 2017), which can be thought of a combination of structural equation model (SEM) (Peters et al. 2017) and vector autoregressive model (VAR) (Sims) 1980). To estimate the SVARM parameters, a stream of rescarch assumes the exogenous noise to be non-normally distributed (Hyvärinen et al.) 2010 Moneta et al. 2013). This allows us lo apply independend component analysis (ICA) to infer the causal structure from observations (Hyvarinen 1999). Then, leveraging non-convex optimization, DYNOTEARS showed promising results the task of causal structure learning for time scries (Pamfil et al., 2020 From the optimization point of view, the alternating direction method of multipliers (ADMM) (Boyd] et al. 2011] is also exploited in several works for causal structural learning, see, e-g-; (Ng & Zhang. 2022:
Yang et al. 2022 Harada & Fujisawal 2021). As detailed in Section 3 also the method proposed in this paper will hinge on ADMM but, differently from Ng & Zhang (2022) and Yang et al. (2022). we leverage a modified ADMM algorithm that exploits a linearization of the non-convex dagness function introduced by Zheng et al. (2018). In Harada & Fujisawa (2021) ADMM applies to a linearly-constrained problem having a different objective function. Also, differently from our method, the approach in Harada & Fujisawa (2021)
is viable only if the latent noise is assumed to be non- non-Gaussian.
Regardless of the assumed causal model, all previous works suffer major limitation, that is, they only assume generative processes where instantaneous and lagged interactions occur at the time resolution corre- sponding to the observational task. However, the introduction of multiscale causal model able to localize causal interaction in both time and frequency, and to provide information regarding the underlying causal mechanism driving the interactions, is of paramount importance since the multiscale property represents a key feature of complex systems, as shown in [Cllet & Fisher (2001); Kwapień & Droždž (2012). Consider for example fMRI data regarding different brain regions of interest (ROIs) collected at timestamps t At, with At being the time interval between two consccutive samples. The aforcmentioned limitation translates into the inability to study the system at hand by means of causal model in which the causal mechanism determining the interactions between ROIs are allowed to occur over time resolutions coarser than At, while it is known that the connectivity between ROIs varies over different time scales, also depending on the state of the brain (Jacobs et al.] 2007 Ciuciu et al.) 2014 Ide et al. 2017). Other examples can be found in other application domains (Besserve et al., 2010: Gong et al. 2015 Runge et al., 2019a), and in general there is no prior knowledge about the time scales at which important causal relations occur.
Noteworthy, researchers active in the neuroscience field have proposed different methods to analyze the functional connectivity between brain areas from causal perspective in the frequency domain [Kaminski & Blinowska, 1991 Guo et al. 2008). All these methods are variations of the well-known Granger causality Granger 1969), and relate to its frequency domain formulation obtained from the spectral decompasition of the VAR process (Geweke 1982 1984). While these methods offer valuable insights into potential causal relationships between signals from spectral analysis, they cannot be considered causal models and do not provide a description of the mechanism of the underlying causal processes driving the interactions, as our model does. The same holds for alternative approaches developed to assess directional connectivity wet
2

---

Published in Transactions on Machine Learning Research (10/2023)
brain signals in the frequency domain by analyzing the distribution of phase differences across different frequency components (Nolte et al. 2008 Stam & van Straaten 2012).
The interest in a multiscale analysis of dependences among time series emerges also from other application domains. In particular, wavelet analysis, together with network analysis, has been already employed in the study of financial risk contagion (Loh) 2013; Khalfaoui et al. 2015 Wang al.) 2017). More recently, the integration of machine learning methods and multiscale representations has been proven to provide significant advancements in biological and behavioral sciences, reported in |Alber et al. (2019) and Peng et al. 2021). Recently, D' D'Acunto et al. (0022) proposed a probabilistic generative model and inference method for multiscale DAGs. In contrast to that work, our methodology (i) relies on wavelet transform, generalizing the SVAR model to the time-frequency domain, (ii) enables learning of multiscale relationships across different time lags while making no assumptions on the causal ordering at each scale, and (ii) has a lower computational cost, thus allowing analysis of networks with a larger number of nodes. For these reasons, we consider the method in D' D'Acutto et al. (2022) incomparable to our algorithm. Overall, the combination of multiscale representation with causal learning is stiill an open problem and this motivates the proposed work.
Contributions. In this paper, we overcome the limitation of previous approaches by proposing a linear causal inference algorithm based on a multiscale representation, able to capture the most relevant causal dependencies across multiple time scales and time lags. Specifically, we start formulating problem aimed at learning a multiscale causal graph from data, while taking into account sparsity and enforcing causality through directed and acyclic topology. Then, as reported in Section 3 we exploit a customized version of lincarized ADMM (useful to copy with the non-convexity of the problem at hand), thus deriving the proposed Mulliscale- Causal Struclure Learning (MS-CASTLF)) algorithm. As a particular case, LIS-CASTLE includes a single time scale version, which we term Single- Causal Structure Learning (SS-CASTLE), to learn causal connections at the frequency of observed data.
optimization
To summarize, the paper has three main contributions:
Multiscale structure learning. Firstly, propose multiscale causal inference algorithm that allows the extraction of causal links at different time scales, without requiring any prior knowledge of the scale where causal relations are most effective. We evaluate the performance of MS- CASTLE from synthetic data generated according to a multiscale causal structure. The results of our empirical assessment, illustrated in Section 4.3 suggest that MS-CASTLE is suitable for both Gaussian and non-Gaussia settings, and that is robust to different choices of the wavelet family used to decompose the input time series.
Application to financial markets. Secondly, we showcase the application of our proposals on real- world financial time series in Section 5 Specifically, we apply MS-CASTLFE to learn the causal dynamics of risk contagion among 15 global equity markets during covid-19 pandemic, from January 2, 2020 to April 30, 2021, and we compare the resulting graphs with those retrieved by single-scale causal learning methods. Our analysis shows that MS-CASTLEE provides richer information regarding the causal structure of the system that cannot be understood by looking only at the estimated single- scale causal graph. In particular, our findings suggest that: i) causal connections are characterized by positive weights and are denser at mid-term time resolution (scales and 4, i.e., 8-16 and 16- 32 days, respectively); ii) the strongest connections are lagged and they appear at scales and 4; iii) the markets injecting the majority of risk into the network are Brazil, Canada, and Italy. Further discussions concerning the obtained results and the richness of information gained through the multiscale causal analysis are given in Section 5.4 Our analysis of financial time series provides novel results at both methodological as well as application levels with respect to the stream of rescarch known as Econophysics (Manlegna & Stanley 1999). At the methodological level, we propose a multiscale machine learning causal model that, differently from existing work (Biliio et al. 2012), allows us to analyze both instantancous and lagged causal interactions at distinct time scales. At the application level, we apply LIS-CASTLE to learn the causal dynamics of risk contagion among 15 global equity markets during c9vid-19 pandemic, rather than focusing on financial institutions.
3

---

Published in Transactions on Machine Learning Research (10/2023)
3.F
Y
Y
It
x
X
X
W12
V
XL
X
4--1
w
S
X
X
x
X
X
(a)
(b)
Figure 1: Example of SSCG with N = 3 and maximum lag order L = 2: b Example of MSCG with N = =3 D time scales and maximum lag order L =
Single-scale structure learning. Finally, we compare the single-scale version of our algorithm, i.e.,
SS-CASTL.E, with several baselines (Pamfil et al. 2020 Runge 2020: Hyvärinen et al. 2010 Hyvarinen
1999 Shimizu al.) 2011) on synthetic data (Sections 4.1 and [.2]. We consider different settings
study the robustness of the performances along (i) the number of available observations, (ii) the size
of the network, and (iii) the distribution of the exogenous noise used to generate data. Our empirical assessment shows that SS-CASTLE compares favorably with all other single-scale basclines, while also sensibly reducing the computational cost.
Notation. We denote by No the set of natural numbers including zero. We indicate with [X], X E N, the range of the numbers from to X and with [X] that from to X. We denote scalars by lowercase letters 2, vectors by lowercase bold letters x, and matrices with uppercase bold letters X. Finally, we denote with || || the Frobenius norm of matrix, with Tr(-) its trace, and represents the Hadamard product.
2 Problem Formulation
Let consider data set Y RTxN composed by N time scries of length T 2D, for some D N. Let yi[4] be the value assumed by the i-th time series at lime 4, and let Pit denole the set of parents (in the DAG representation) of yt[t), with lag I E No. In the single-scale causal structure learning problem for time series, we are interested in understanding whether the considcred time scries admits a functional representation in which yr[ depends on a set of parent variables, up to a finite lag L:
y: [t] = f (P....... PP......PPPP.....
i [N],
(1)
where c(t) represents either additive noise, statistically independent of the i-th time series, possible model mismatch, occurring at time t. It is worth noticing that the set of parents Pt, vary with t, To distinguish causes from effects, the system of Equations ( must admit representation based on a DAG. As far as lagged interactions are concerned, since we cannot observe causal elfects from present to past, P., may contain (t - 1J, with / > 0. In other words, time ordering provides lagged causal connections with implicit causal direction. However, when we look at instantancous interactions, if we represent each P., over graph, then the overall graph must be acyelic, otherwise it would be impossible define the direction of the causal relation.
4

---

Published in Transactions on Machine Learning Research (10/2023)
By limiting our attention to linear dependencies, the causal inference model can be expressed as:
y|t] Ext ylt lW e(t]]
(2)
which coincides with the so-called SVARM. In Equation (2) y(t] := "uldf.....,y [t) E RIxN is the row vector containing the values assumed by N time series, at time t; whereas W E RNXN, with I E [4]o and being the maximum lag, is the matrix representing the cansal relation at lag 4, so that w/y * 0 if yi(t - -] E Pa.. In particular, W represents instantaneous interactions and its structure is such that, if we map the coefficients of W over the edges of a graph of size N, the resulting graph is acyclic. Finally, <|1] E RIXN represents a random disturbance or model mismatch al time Equation 2 is said to be structural since it allows us to express variables (effects) as linear functions of other exogenous variables (causes), considering instantaneous as well as lagged relations, also referred to as intra- and inter-layer connections, respectively. As example, Figure Ia shows the single-scale causal graph (SSCG) associated with Equation ( in case of N = 3 and T= 2. In Figure Ta the subscript represents the node index while the time lag is given within the square brackets. Notably, the graph represents causal interactions from the past to the present, with instantaneous effects at time ( avoiding cycles.
2.1 Multiscale Causal Inference Model Based on Stationary Wavelet Transform
The model represented in Equation ( and sketched in Figure la although very well known and applied, is however limited because it implicitly assumnes that the time scale at which important dependencies show up is that associated with the obsorvation task. As discussed in Scction [ cansal dependencics might occur at different time scales. Hence, an intriguing and underexplored research problem involves the scenario where the observed variable ys[(] in Equation ( can be represented by a multiscale functional decomposition. This means that y: may be decomposed into contributions vf [t), with d E [D], each having a representation in the form described in Equation (1) and related to certain frequency band. Mathematically, by defining P as the parent set of (1(4) at the d-th time scale, we can express y (i((] as
D 3:14 = Σ gid (Pr......, Pr Pi... e.al()) d=1
(3)
where (ia[(]] represents an exogenous noise associated with the d-th time scale. Thus, in the case of linear multiscale causal dependencies, it would be beneficial to enhance the SVARM with a multiscale modeling approach.
The first question to tackle how to obtain useful multiscale representations for modeling causal relationships that exist the time and frequency domains. To this aim, among the various existing data transformations, we have chosen the wavelet transform. There are several reasons underlying this choice. First and foremost, this transformation, when equipped with orthogonal filter family, allows the decomposition of input signals without loss of information. This is a key feature that enables a perfect reconstruction of the signal, in accordance with Equation (3). This property is closely related to the energy preservation guaranteed by this transform, which preserves the original interpretation of the signal as it does not amplify or suppress different characteristics of the signal's power.
As a sccond motivation, the wavelet transform allows the localization of a signal in both the time and frequency domains. This is because the filters used by this transform have compact support, meaning they are nonzero only over finite interval. Therefore, the values of the obtained representations will depend on the information contained in respective time intervals determined by the filter used. The time- frequency localization, which distinguishes the wavelet from the Fourier transform, is essential for our work as Equation assumes interactions that might occur at different time scales and with different time lags.
Lastly, the wavelet transform is non-parametric transformation. This final characteristic makes it applicable to wide category of signals, as the transform does not assume any specific functional form of the signal to be processed, nor does it rely on particular distributional assumptions.
5

---

Published in Transactions on Machine Learning Research (10/2023)
Hence, to develop our multiscale causal inference model, we apply wavelet decomposition to the observed data set Y. In the sequel, we provide only the essential information about the framework; more details are given in Appendix A] Each row y[t] of the data set represents a sample collected at the timestamp t. At, where At is the time interval between consecutive samples. The wavelet decomposition of level transforms each time series Yt into - vectors of wavelet coefficients and an additional vector of scaling coefficients. The d-th wavelet coefficients vector corresponds to the variations of y; at time scale 26-1, At, representing the frequency band [1/24+1, 1/24]. These wavelet coefficient vectors capture the input signal variations over time scales ranging from At to 2D-2 At, corresponding to frequencies from 1/2" to 1/2. The scaling coefficients vector contains information about variations over the scale 2D and coarser scales, representing frequencies slower than 1/2D. Our notation assigns the finest scale (d = 1) and the coarsest scale (d (d= D) accordingly.
While the proposed model flexible and can accommodate various types of wavelet decomposition (Percival & Mofjeld 1997), our study found that the stationary wavelet transform (SWT) (Nason & Silverman 1995) has the best. performance for the learning task. In our notation, [ E RN represents the SWT at time and scale d for the N time series, where d E [D] and D corresponds to the scaling coefficient. The SWT provides non-decimated detail coefficients [t) at each scale, offering a translation invariant representation. This property is advantageous as it captures relevant information without considering the position of the analysis time window. In addition to the motivations above, to preserve both odd and even decimations at each decomposition level and avoid unnecessary redundancies, we consider orthogonal filter families within the SWT framework. This choice ensures that the time series yf associated with different time scales are orthogonal. Furthermore, due to the orthogonality property and as previously discussed, the energy of the input signal is conserved by the transform and distributed across the scales (Percival & Walden 2000).
Collecting D different scale details into single vector y := y"-'(, .fl.....y'l y [t]] and stacking these row vectors for t [T] on side of cach other, we build an augmented data set Y E RTxN with N = DN. In this matrix, the row vector at timestamp L contains the 1-th detail values of scale d for all N input d + 1) N. Then, we build the block diagonal matrix W block[WP, W ....., Wi] of size Ñ x Ñ. Each d-th block W in W represents the causal interactions d. occurring at scale d with lag I E [T]u, where w/y * 0 if and only if vf - 1] E Py, Here, Py, represents the parent set for time series y; at lag / and scale d, which can vary across both graph layers and pages. Then, the resulting multiscale causal inference model reads as:
signals, indexed as (D d) N + +1, 1......
wij
y = Σ yl - W + Z[0 1=0
(4)
where 2[t] E RIIN denotes the additive noise term. In the proposed model, the matrix wo must satisfy the acyclicity requirement lo set causes apart from effects.
A pictorial example is given in Figure lb which depicts a multiscale causal graph (MSCG) for N = 3 time scries, D = 3 time scales, and maximum lag L = 2. Each layer in the graph corresponds to a specific time lag, while different pages represent different time scales. In the notation used in Figure 1b the node superscript denotes the scale index (page of the graph), and the subscript indicates the node index. The time lag is indicated within square brackets. Within each page, we observe both inter-layer and intra-layer directed edges, where the latter, similar to the SSCG case, do not form cycles. However, variables may exhibit different interactions at each time resolution. Therefore, when considering interactions across pages, reverse causal relationships between variables can be observed, as shown by the blue and orange ares in Figure [b] Additionally, due to the use of an orthogonal wavelet family, there are no ares between pages. Thus, Figure [b] represents a multiscale DAG that incorporates both instantaneous and lagged linear causal relations at different time scales. Each page of the graph is an SSCG at a specific time resolution.
2.2 Optimization Problem Statement
The multiscale causal inference problem aims at learning the matrices W E [4). in Equation 4 from data. To ensure the acyclicity of the estimated MSCG, the inferred matrices of causal effects W must
6

---

Published in Transactions on Machine Learning Research (10/2023)
entail DAG. However, learning DAGs from observational data is combinatorial problem and, without any restrictive assumption, it has been shown to be NP-hard (Chickering et al. 2004). In our case, since we cannot observe edges coming from the present to the past, lagged causal relationships encompassed in the matrices W (with > 0) are acyclic by definition. Therefore, the main issue concerns the inference of the matrix W representing instantaneous causal effects. To handle the acyclicity of W similarly to DYNOTEARS Pamfil ct al. 20200,, we exploit the dagness matrix function h(M) VNN + R proposed by Zheng al. (2018). who proved that matrix W of size N N can be represented DAG if and only if
h( k(W9) = Ir
-N
N
(5)
Using this function as penalty term, non-convex, optimization problem. To this aim, let us introduce the matrix Y [Y Y Y4 RTV (with V N(L + 1)) containing the matrices of I-shifted observations Y Rxx Similarly, we build W = W W WT" RVXN, For convenience, let us indicate with B the set of matrices having the same structure as W, ie., made up of stacked block diagonal matrices. Then, the proposed multiscale causal structure learning problem can be mathematically cast as:
able formulate the learning task as continuous, albeit
min WEB
f(W) + W||
(6)
subject to h h(W") Tr (e (eW%owoo
- === =0
where
/W)= = 1 - YW
The & norm penalty used Problem (6) enforee sparsity of the aggregated causal matrix W, with a tunable parameter À > 0. Here, differently from already proposed ICA-based estimation procedures (Hyväii- nen et al. 2010: Moneta et al. 2013), the matrices of causal coefficients learned simultaneously. Despite the convexity of the objective function, Problem 6 is non-convex due to the presence of the acyclicity
constraint h( (W) 0. In the next section, we derive an efficient method to solve Problem (6)
3 The MS-CASTLE Algorithm
To find local solution of Problem (6) or, more precisely, point that satisfies the Karush-Kuhn Kuhn-Tucker (KKT) conditions, we exploit the computational efficiency of ADMM (Boyd et al. 20011). To this aim, we recast. Problem 6 the following equivalent manner, introducing the auxiliary matrix Z RVXN and obtaining
min WeB
f(W) + 1
subject to h (W") Tr (T (wWowwo N
(7)
w -Z=0vxx
Let us now denote by W = vee W) E RRV and z = vec(Z) (Z) E RVR Then, following the scaled ADMM approach (Boyd et al.) 2011). and letting a and & be the Lagrange multipliers associated with the equality constraints of Problem we introduce the following augmented Lagrangian (AUL) function:
Cp (W2 f(W) + ah w A1|2/11 + l ( 2.
(8)
where g(W,,, 2, = 2 - 2 + 013 2 1111. and p 0 is a tunable positive coefficient. The ADMM
2
algorithm proceeds by iteratively minimizing the AUL function with respect the primal variables W, z, while maximizing it with respect to the dual variables a and B, looking for saddle points of the AUL. However, while the AUL function is strongly convex w.r.t. 2, and naturally concave w.r.t. and B, it is non-convex w.r.l. W, due to the presence of the non-convex dagness function (Wo) in Equation (8) To handle this
7

---

Published in Transactions on Machine Learning Research (10/2023)
non-convexity issue, following the idea of linearized ADMM methods (Yang & Yuan 2013 Goldfarb et al. 2013), we substitute the non-convex dagness function h (wo) in the AUL with its lincarization around the current value W assumed at each iteration k, i.e.,
h (W" w?)
(w)) + IT (G(W)) w!))
where G(W") represents the matrix- gradient of function h (w"). Then, substituting Equation ( we obtain the following approximated AUL:
9 into
P Ip W. (W,,,,,::W)) f(W) + ah (W:WW)) + A||z|11 2 2, 3),
(10)
which is now strongly convex w.r.t. W while preserving the first-order optimality conditions of the AUL in Equation (8) around the current approxination point W As a result, any point satisfying the KKT conditions using the approximaled AUL in Equation (10). satisfies also the KKT conditions of the original Problem (7). Hinging on this fact, we now apply ADMM to the approximated AUL in Equation d10p. Then,
letting W k. ZE, Ok, and B& be the current guesses of the primal and dual variables at time k we obtain the following set of recursions:
W K++1 arg min f(W) + on Ir (G(W2) Wo + -. Zu + Bx
(11a)
WeB
Zk+1 = arg min Al|z|| + W 2 3x
(11b)
0k+1 = 0 + + h W k+1
(11c)
8x+1 = Bg + Wk+1 Zk+1
(11d)
The first step in Procedure is the minimization of a strongly convex quadratic function, subject structure constraints W B, i.e. simple linear constraints on the elements of W. We perform this mini- mization using the L-BFGS-B algorithm (Byrd et al. 1995). i.c., variation of the Limited-memory Broy- Sen-Fleeeher Golddarb-Solaa Shanno method that handles box constraints. The second slep in Procedure
can instead be computed in closed form as Boyd et al. 2011
2b+1 S(A/p) (W++1++ + B
(12)
where S5(x)= = sign(x) max(# 0) the element-wise soft-thresholding function, used to enforce sparsity of the causal matrix representations. The third step in Procedure (11) performs gradient ascent step to maximize the function in Equation 110) with respect to a, using a (possibly line-varying) step-size 7.
Similar arguments then hold for the fourth step of Procedure 11 All the steps are then summarized in Algorithm 1, which we term MS-CASTLE.
Remark:
Algorithm D can be easily customized to solve Equation 2 ,, where we simply ignore the multires-
olution analysis. With regards Algorithm 1, simply means skip line 2. This leads to the aforemen- tioned SS-CASTLE algorithm, which applies a particular sub-case of Problem 6), in which we have: (i) WL E (RN(++))N,, where
Y Y; (ii) Y = (YO, Y Y!...,, Y E Tx+(x(++ (iii) W = wwo" WI W' E NNN are the matrices of causal coefficients of Equation
Regarding the initialization of W and Z, they must be initialized in order to satisfy the primal feasibility conditions of Problem 6 i.e., to satisfy the two constraints. A possible choice is to set them equal to zero. Regarding the dual variables of the augmented Lagrangian in Equation ( and B, since they concern equality constraints, they do not need to satisfy any specific condition. Thus, they might be initialized to zero as well.
Initialization.
AUL parameters. The papee--araneters p and y represent the augmented Lagrangian penalty parameters that be viewed as the equivalent of dual-ascent. step sizes in the ADMM procedure. Larger values of
8

---

Published in Transactions on Machine Learning Research (10/2023)
Algorithm MS-CASTLE
procedure MS-CASTLE(Y, L, p,9,,,,,," maxite)
Y Apply SWT Y
Initialize ,...,
Y + I Y- Y4
5: while k < maxiter & h4 > do
6:
Find WK+1 solving Problem 1la
7:
18+1 + -h (W(-))
8:
if hx+1/hg > r then
Y + 7
D. Y E (0, b Soft-threshodding
10:
Zk+1 + S(A/p) (Wi+1 + (x)
11:
0k+1 + 0g 7 h&+1
12:
B t Bx + Wk+1 Zk++
13:
return W
these step sizes result in greater penalties for constraints' violations (i.e., primal feasibility). In our case, in (Zenng et al. (2018), the increase in y is practical way to manage the violation of the acyclicity constraint. Consequently, it determines the speed of convergence of the proposed algorithm.
Computational cost per iteration. Let us consider w.l.o.g. the single-scale case. In fact, by considering only the active parameters to optimize, the cost of the multiscale case is simply D times that of the single- scale. To solve Equation (1la), we apply L-BFGS-B B to an objective function that has cost of O(N2) operations when T 2 N and L < N. In fact: (i) S(W) requires the computation of the product YW (which costs O(N2)), plus TN subtractions and the calculation of the squared Frobenius norm (which costs O(TN)); (ii) Tr WW0 requires O(N2) operations to calculate the product and O(N2) operations for the trace; (iii) ||w zk + Akll requires O(N2) additions and subtractions plus O(N2) operations for the squared fz-norm. To solve Equation (IIb,, we apply the element-wise soft-thresholding operator vector of size (L 1). Therefore, the cost is O(N2) operations. To solve Equation {11c). need to compute the dagness function of an N X N matrix, which has cost of O(N3) Zheng al. 2018 Finally, solving Equation 11db costs O(N2) operations.
In the next section, we illustrate the performance of MS-CASTLE and SS-CASTLE, comparing them with alternative methods available the literature.
4 Numerical Results
In this section, we start showing the advantages of SS-CASTLE (i.e., the customization of the proposed MS-CAS CASTLA method to temporal causal structure analysis) over existing alternative methods solving Equation 2 More specifically, Section 4.1 shows that, when compared DYNOTEARS (Pamfil et al., 2020), which aims to solve the same optimization problem, SS-CASTLE benefits from the linearization procedure described above to reduce the computational cost of each iteration while preserving performance. In addition, Section 4.2 illustrates how SS-CASTLE outperforms non aan-gradient--bee
methods, e.g., VAR- ICALiNGAM, VAR VAR-DirccLLINGAM (Hyvärinen et al. 2010). and PCMCI+ (Runge) 2020), when we sample <[(] from a p-generalized normal distribution, with p E {1, 1.5, 2, 2.5, 100}. Here, we evaluate the behavior of SS-CASTLE along the size the network N and the number of samples as well. Finally, in Section 4.3 we assoss the performance of MS-CASTLE, on multiscale synthetic data. Also in this case, we test the performance of our method along N, T, and both Gaussian and non- non-Gaussian settings. In addition, we study how different choices of (i) wavelets and (ii) X hyper-parameter affect the performance of MS-CASTL.E.

---

Published in Transactions on Machine Learning Research (10/2023)
Number of iterations
Time (9), log scale
10a
DNROTEARS SS-CAbILe
104
1. 12 12.
10
10
50
100
30
100
Nodes
Nodes
Figure 2: Number of iterations (left) and computational time (right) to solve Equation shown by DYNOTEARS (purple) and SS-CASTLE (pink).
4.1
Comparison with DYNOTEARS
When reduced to the single-scale temporal version, Problem G corresponds to Problem (5) in Pamfil et al. 2020 where the l regularization paramcters for lagged and instantancous causal relations are set to be equal. However, unlike the DYNOTEARS algorithm proposed by Pamfil al. (2020), which solves the optimization problem using augmented Lagrangian and relies on the solution scheme proposed by Zheng et al. (2018) to handle non-smooth terms, SS-C BS-CASTLE provides solution by leveraging linearized ADMM. Specifically, from theoretical point of view, considering the analysis of computational cost per iteration of our method reported in Section 3 we have that SS-CASTLE reduces the compulational cost of DYNOTEARS in two ways. Firstly, the linearization of the dagness function lowers the computational cost required to update W. In fact, in our case, the calculation of the objective function requires O(N²) operations, compared to the O(N3) required by DYNOTEARS. Secondly, to impose regularization on the & norm of causal coefficients, CSSCAE leverages the element-wise woft-thresholdig operator. In contrast, similar to tooZengg et al. (2018), DYNOTEARS uses the splitting trick, writing W W+ W-, -, where W+ 2 0 and W_ 2 2 (having the same dimensions W), thus doubling the number of parameters to be estimated. In the following sections, we further expand the comparison between SS-CA CS-CATLLE and DYNOTEARS. We test the two methods over synthetic data, so that the ground truth is known. Our goal is to compare empirically the computational time needed to the two alternative methods to estimate the causal matrices in Equation with similar accuracy.
(2)
4.1.1 Data Generation
We generate synthetic data according to Equation 2 Specifically, we set I = and we assume that each ei[t] ~ N (0,o?) with of E [1, 2]. Moreover, we set the number of samples T = 1000 and we pick NE {10, 30, 50, 100). For each of the four possible values of the number nodes, we simulated 100 data sets. Regarding the causal matrices, we generate them by adopting the same procedure illustrated by Hyvärinen et al. 2010). To manage the level of the sparsity of Wo and W we introduce the parameter s E (0, 1). The latter is used as parameter of a Bernoulli distribution, more precisely B(1 s), which controls the number of nonzero coefficients of the causal matrices. As the number of nodes N grows, we increase the sparsity of the causal structure. More specifically, the combinations (N, used in the experiments below are ((10, (30, .85), (50,.90) (100,.95)):
4.1.2 Results
Figure 2 displays the mumber of iterations (left) and computational time (right) needed to solve Equation function of the number of nodes, by DYNOTEARS (purple) and SS-CASTLE (pink), respectively, On the left of Figure 2 we report a swarm plot in which, given a certain number of nodes, each point represents the number of iterat ions required by each algorithm to retrieve the solution. In accordance with the data- generating process described above, for each value of N we have 100 points per algorithm. Therefore, given a certain value of N and a specific number of iterations n, the number of points reported horizontally represents the number of data sets (composed by N time series) in which the algorithm required n iterations solve the problem.
2
10

---

Published in Transactions on Machine Learning Research (10/2023)
On the right of Figure we provide violin plot that depicts, for each value of N, the histogram of the computational time (measured in seconds) needed by each algorithm to solve the problem. Moreover, dashed lines within the histogram represent quartiles. From Figure 2 (left), we see that, even though SS- CASTLE needs more iterations to converge, SS-CASTLE significantly reduccs the overall computational time to converge. Furthermore, we also observe that the larger the network size, the greater the gain. Hence, in accordance with the theoretical insights highlighted above, we conclude that S5-CASTLEE decreases the computational cost associated with each iteration while preserving performance (see Appendix B
4.2 Comparison with Non-gradient based Methods
We compared the performance of SS-CASTLE with two major linear non-Gaussian methods, namely VAR- ICALiNGAM and VAR-DirectLiNGAN
well as a constraint-based method called PCMCI+ synthetic datascts. More in detail, the former two models rely upon the assumption of non-Gaussianity of c(t) in Equation VAR-ICALINGAN belongs to the family of ICA-based methods: first it fits a VAR model to recover lagged causal interactions and then it employs FastICA (Hyvarinen, 1999) on VAR residuals to uncover instantaneous relationships. In the past, several ICA-based algorithms have been developed. How- ever, as shown in Moneta et al. (2020), previous models are equivalent in terms of performance. Regarding VAR-DireetLiNGAMM it was proposed in order to solve the possible convergence issues of ICA-based meth- ods (Himberg et al., 2004 and it is guaranteed to retrieve the right solution of the problem if the model assumptions are satisfied and the sample size is very large. Concerning PCMCI+, it represents an improved version of PCMC1 (Runge et al.) 2019b), and is able to deal with both linear gngear/nonle nonlineer lagged and instanta- ncous causal relationships under the common assumption of causal sufficiency, faithfulness, and the Markov condition (Peters et al. 2017). In the experiments below, in order to the aforementioned models, we use the lingam and Tigramite Python package=" made available by the authors.
4.2.1 Data Generation
We generated synthetic data by using Equation P in which we set L 1. Moreover, we conducted an extensive simulation study order to assess the robustness of all the methods in different. settings. Tn particular, we varied the features of the generated data sets as follows. Firstly, we use different data set sizes, T E {100, 500, 10003. By varying the number of samples, we can inspect the sensitivity with respect to the data set size of the tested algorithms. The latter aspect is relevant in several fields, especially when the system at hand shows non-stationarity. For instance, this is the case in finance, where practitioners usually deal with a small number of historical observations due to the continuous evolution of financial markets. Secondly, we vary the network size, N E {10, 30, 50}. Concerning the level of sparsity and the generation of the causal matrices, we adopt the same methodology described in Section 4.1.1 Last but not least, we sample E[t] from p-generalized normal distribution, with p E {1, 1.5, 2, 2.5, 1003. The p-generalized normal distribution is defined as follows (Kakke & Richter 2013).
Definition 4.1 (p-generalized normal distribution). Let us consider x E R, P R+ Therefore, the p- generalized normal distribution has density function equal to
p 1/1//D 1/D
fz(z)
21 (1/p)
exp
P
(13)
where I is the gamma function.
The parameter determines the rate of decay of Equation
More in detail: (i) 1 corresponds to a [ii for p = 2 we get the normal distribution; (iv)
13
Laplace distribution; (ii) p= 1.5 the super-Gaussian case; p == 2.5 is the sub-Gaussian case; (v) for p = 100 we obtain approximately a uniform distribution. Therefore, as p diverges from 2, the non-normality of e(t) is enhanced. For each combination of the parameters, we generate 100 data sets.
11

---

Published in Transactions on Machine Learning Research (10/2023)
7 100
T=
1000
N
N =0
1.01
1.0
022
0.0
1.01
1.0
PMMM
PONEL 55 CASTLE
95. CASTES
0.04
1.0
1.5
2.0
25
10h.
10
1.5
20
25
102.0
1.0
1.5
20
25
1000
100
1.5
2.0
25
105.0
1.0
15
20
25
100.0
1.0
1.5
20
25
100.
(a)
(b)
Figure 3: Comparison with the baseline methods in the estimation of the causal matrices w (top) and W
(bottom). Each subfigure depicts box plots for the F1-score when a (N, T) [30} (100, 500, 1000) and b (N, T) E x {1000}.
4.2.2 Results
Before testing SS-CASTLE on the generated data, fine-tune the sparsity strength À onto separate data sets generated according to the procedure explained above. More precisely, we let À assume values in the set (000,,000,, 00,, and, for each combination (I', N,,N, we chose the best value according to Fl-score and SHD (structural Hamming distance). Please, refer to [Appendix C for the definition of these metrics. Specifically, we first compute two rankings of values. In the first ranking, we sort the values of À according to the Fi-score descending order (the higher, the better). In the second, we sort the values according to SHD in ascending order (the lower, the better). At this point, for each À we sum the positions obtained in the two rankings, and we pick the X showing the lowest value. Due to the needed computational time, for the case N = 50 restrict the possible values of X to {.01, 1, .5.. The latter restriction does not impact the analysis since our objective is not to find the optimal value of the hyper-paraneter, but rather to set the latter in a dala-driven manner. The chosen values for À are given in Appendix D
Sensitivity to data set size. Figure 3a depicts the performance comparison in the learning of Wo (lop) and W (bottom). For readability, we provide only the results in the case of N 30 nodes. Results are qualitatively equivalent in the other two cases. As can from Figure Bal SS-CASTLE outperforms VAR-DireeLiNGAM, VAR-ICALINGAM, and PCMCI+ in terms of Fl-score. In addition, Appendix E provides results for SHD metric. Concerning the learning of W even though VAR-DireeT.iNGAM shows a slightly better performance than VAR-ICALINGAM in case of strongly non-Gaussian settings (p = 1 and p = 100) and larger data sets (T = 500 and T = 1000), it tends to suffer more when the Gaussaasiyy assumption becomes violated and a lower number of samples is available. The latter behavior is consistent with DirectLiNGAM model assumptions. Regarding PCMCI+, its performance does not vary along p. consistently with its formulation. However, it tends to underperform the rest of the tested methods. In accordance with the problem formulation (see Section 3), SS-CASTLE does not show any dependence on p, and requires a smaller number of data to converge towards more accurate solution.
Regarding the matrix of lagged causal effects W differently from the considered non-Gaussian methods, SS-CASTLE contemporaneously estimates inter and intra-layer connections. In addition, we see that all models tend to be more accurate in retrieving the lagged interactions. Please notice that, given the time ordering, W is acyclic by definition. Therefore, all entries could be different from zero. This means that in the case of N = 30 and 0.85, on average we have 135 nonzero coefficients.
Sensitivity to network size. Figure 3b shows the comparison of models performance in the estimation of the causal relationships along N. For readability, we provide only the results for the case T' = 1000. The results are qualitatively equivalent in the remaining two cases. In addition, Appendix provides results for
The packages are available
and http::///sshaabcccc////aannngggr
12

---

Published in Transactions on Machine Learning Research (10/2023)
SHD metric. Looking at the case of instantaneous interactions (top row), at first glance we notice the higher error variance in the box plots related to F1-score, when N = 10, across all models. However, even though the Fl-score is a normalized metric, this is an effect of the small number of instantaneous causal connections. Indeed, in the case of N = 10 and & = -80, on average the ground truth W has only 9 entries different from zero. As a consequence, a single mistake weighs more. Overall, from Figure 3b we can see how SS-CASTLE outperforms the other methods. Moreover, we do not appreciate a decrease in performance when the number of nodes increases. In addition, SS-CASTLE proves to be robust changes in the value of p.
Regarding lagged causal connections (bottom row), the models tend to perform better than in the case of instantancous relations. We underline that, in the case of N = 10, we do not observe the same error variance as above. Indeed, with the same level of sparsity 8, on average W has approximately twice as many non-zero entries than w Overall, we notice that baseline algorithms tend to suffer when the mumber of nodes increases. Furthermore, our results show that VAR-ICALINGAM tends to be more robust than VAR- DiMM as the non-Gaussianity assumption turns out to be violated. Also in this case, SS-CASTLE does not show any decrease in performance while varying p. In addition, il achieves high performance in the case of larger networks well.
4.3 Performance of MS-CASTLE on Synthetic Data
Here we provide the empirical assessment of MS-CASTLEE. Analogously to the previous section, we evaluate the performance of our method when varying the network size, the data set size, and the distribution of the random disturbances. In addition, we are interested in understanding how the choice of a wavelet in the decomposition phase affects the performance of MS- CASTLE Therefore, we compare different versions of MS-CASTLE. First, we consider the case in which MS-CASTLEE receives as input the perfect decomposition of the time series, i.e., Y(4 used in the data generation process. This case represents an upper bound for our method, as there are no errors induced by the wavelet transform. We refer to this version as MS- CASTLE opt. dec.. Next, we consider the case in which MS-CASTLE uses the same wavelet used in the generative process, namely Daubechies with a filter of length 10 (MS-CASTLE sut db5). We then test two other versions in which we use a different wavelet (but belonging to the same family) than the one used in the generative process, namely Daubechies with filters of length 6 (MS- CIS-CASTLE sut db3) and 14 MAS-CATTLE sut db7). Finally, we include two additional orthogonal wavelets belonging to different families, the symlet with a filter of length 10 (MS-CASTLE swt sym5) and the coiflet with a filter of length 12 CMSCASTLE sut coif2).
Lastly, we examine the impact of the hyper-paramett choice on MS-CASTLE performance. Additionally, we present a methodological approach tailored for scenarios where additional data sets are unavailable for fine-tuning X.
4.3.1 Data Generation
We generate multiscale observable data Y E RTxN determined by an underlying multiscale causal structure according to the proposed model in Equation 33 and Equation 44. First, we sample an i.i.d. noise e[t] E RN from a p-generalized normal distribution. Next, w.l.o.g. we decompose <[t] with a DWT using as wavelet Daubochics with filter length 10 and considcring = 4, thus obtaining č[t in Equation (4). Hence, we generate y by imposing a multiscale causal structure containing both lagged (L= 1) and instantaneous interactions four different frequency bands (see Appendix F). Finally, we obtain y/t] by summing the contributions from the time scales, in accordance with Equation (3) and in general with the synthesis capability associated with the wavelet transform (Percival & Mofjeld 1997).
In the experiments, we perform the following analysis: (i) we study the sensitivity to the size of the graph we vary N E (5, 30); (ii) we assess possible dependencies on the number of samples we vary T e {128, 512, 8192]; finally (iii) we deal with the Gaussian p = 2 and non- non-Gaussian p = 1 cases. For each of these settings, we simulate 20 data sets.
13

---

Published in Transactions on Machine Learning Research (10/2023)
N 5
N 10
N = 30
1.0
Gaussian case
01.8
0.6
1.0
non- non- Gausians case
0.8-
MSCASTLE <pt. dec
MSCASTLE did
0.6
MSCASTLE owt dbS
MSCASTIF wet dh7
MSCASTLE
MSCASTI owt ool?
128
512 T
8192
128
512
8192
128
512 T
8192
T
Figure 4: Comparison of different versions of MS-CASTLE in terms of F1-score, obtained in the experimental settings described above. Subplots on the top refer to the Gaussian setting, whereas those on the bottom to the non-Gaussian one. Subplots on the left refer lo the case N = 5, those on the center to N = 10, and those on the right to N = 30. Dashed lines represent the quter-quartile range.
4.3.2 Results
Figure depicts the F1-score obtained by the considered versions of MS-CASTLE on the settings described above. The values of A used in the experiments are given in Appendix D while Appendix Grroviies other results concerning the SHD. Overall, MS-CASTLE achieves the best performance when reccives as input y[1] used in the generative process (perfect decomposition of the observed time series). Furthermore, for the same network size N, the Fl score improves as the number of samples T increases. In contrast, we do not notice any dependence of MS-CASTLE on either the number of time series N or the noise distribution (Gaussian vs. non-Gaussian setting). In addition, the results show that the performance of MS-CASTLE is robust to the choice of wavelet used to decompose the input time series.
Sensitivity to À hyper-parameter. Figure 5 shows the behaviour of the F1-score along the regularization parameter X, obtained in the experimental settings described above. For readability reasons, here we show only the results concerning MS-CASTLA swt db5. The plots concerning the rest of the tested models are given in Appendix G. Results show that in case T {128, 512}, adding a sparsity regularization associated with À >> 10-4 effectively improves the F1-score of the model. Conversely, if T is large, there is no benefit in increasing the regularization strength. Obviously, when the value of X becomes too large, the regularization term in the objective function of Problem (6) dominates the first term, i.e., the fitting loss, leading to worse performance.
Network of highly persistent edges. When dealing with synthetic data, it is possible to fine-tune the value of À due to additional data set availability and knowledge of the ground truth. However, real-world scenarios this is not possible, hence a procedure to choose the value of À is needed. In particular, due to the non-convexity of Problem 6 MM-CASTLE can generally converge to stationary points that, possibly, can be different from each other for diverse values of the pparsity-inducing parameter X. Thus, to reduce this ambiguity, in our analysis we look for solutions of the MSCG that are as persistent as possible with respect to different values of A. To this aim, we first choose a suitable range for the previous hyper-parameler, looking at the regularization to fitting loss ratio, i... the quotient of the division between the second and the first term of the objective function of Problem (6). Once the range of is identified, we define the persistence of causal relation a threshold è as follows. Let us indicate with the vector constituted by the regularization to fitting ratios r% corresponding to the chosen k values of X. In addition, consider wi
14

---

Published in Transactions on Machine Learning Research (10/2023)
T 128
T 512
T 8192
1.00
Gaussian case
0.75
0.50
0.00
1.00
non Gaussan case
0.75
0.25
0.00
10
10-2
10
10
10
3
10
10 10-2 10
Figure Behaviour of Fi-score along the considered values of X, obtained in the experimental settings described above. Subplots the top refer to the Gaussian setting, whereas those on the bottom to the non- Gaussian one. Subplots the left refer to the case T 128, those the center to T 512, and those the right to T = 8192. Dashed lines represent the inter-quartile range, while colors refer to different values of N.
the causal coefficient from node to estimated for À = de. Then, the persistence of the causal coefficient is
Σe 1u0 Maotiaa TR
Pij
Σs Tk
(14)
where Juiiis equal to iff |uty] > c, and zero otherwise. Equation 14 assigns higher persistence value arcs that are present in causal structures estimated from Problem for multiple values of X. Also, from Equation (14b, it holds Pi E [0, 1]. However, the formula does not provide any guarantee regarding the stability of the sign of the causal relation. Indood, it only considers the presence of an are and not the value (and therefore the sigu) of the causal coefficient associated with the arc. Thus, here we define as highly persistent only those edges, with Pij > 0.95, which show a stable sign of the corresponding causal coefficient for all values of A.
According to the procedure above, must choose proper range for X. Let consider the non-Gaussian setting above with (N, 7) = (10, 128), which is also relevant for the financial case study discussed in Sec- tion 5 Figure 0 shows the behavior of regularization to fitting loss (on a logarithmic scale for readability reasons) along the considered values of X, obtained in the experimental settings described above. Look- ing at Figure we see that the more performing values of are 0.05 and 0.1. Tn particular, as depicted in Figure these values fall within the range [0.1, of the regularization to fitting loss ratio. Hence, we consider values of associated with regularization to fitting ratio falling into the latter range, i.e., X {0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1, 0.15). Then, retrieve the network of highly persistent edges using 2 = 0.1.
depicts the comparison in terms of F1-score and false positive rate (FPR) between the MSCG obtained by using the best À among the considered, and the MSCG of highly persistent edges. Looking at the boxplot on the left, we see that the network of highly persistent odges underperforms in terms of F1-score. Nevertheless, the boxplot on the right shows that the latter network leads a much greater reduction of the FPR, thus providing
Figure 7
robust. estimate.
15

---

Published in Transactions on Machine Learning Research (10/2023)
128
T = 512
T 8192
Gaussian case
10-1
10-2
10-3
non- Gaussan case
10-1
10-2
10-3
10
10-3
10
10
10
10-3
10
10
10
10
10-2
10-1
X
Figure 6: Behaviour of regularization to fitting loss (on logarithmic scale) along the considered values of X, obtained in the experimental settings described above. Subplots on the top refer to the Gaussian setting, whereas those on the bottom to the non-Gaussian one. Subplots on the left refer to the case 128, those on the center to = 512, and those on the right to 8192.
F1
FPR
0.70
0.65
0.06
Value
0.60
0.04
0.55
0.50
0.02
0.45
MSCASTLE
MSCASTLE pers.
MSCASTLE
MSCASTLE pers.
Comparison in terms of Fl-score (left) and false positive rate (right) between the MSCG learned by using MS-CASTLE sut dbō provided with the best value of À and the MSCG of highly persistent edges. The results refer to the non-Gait Con-Gaussian setting, where N 10 and T = 128.
Figure 7:
5 Causal Structure Analysis of Financial Markets
In this section, we apply the proposed technique to infer the causal structure of financial markets. We consider data concerning 15 global equity markets at daily frequency. To focus on covid-19 19 pandemic period, we restrict our attention to observations from January 2020, the 2nd to April 2021, the 30th. In our analysis, deal with the following markets: All Ordinaries Inder (AOR, Australia), Hong Seng Indez (HSI, Hong Kong), Nikkei 225 Inder (NKX, Japan), Shanghai Composite Indez (SHC, Shanghai), Straits Times Inder (STI, Singapore), TAIEX Inder (TWSE, Taiwan), DAX Inder (DAX, Germany), FTSE MIB Inder (FMIB, Ital)), IBEX Inder (IBEX, Spain), CAC 40 Index (CAC, France), FTSE 100 Inder (UKX, UK), RTS Inder (RTS, Russia), Bovespa Index (BVP, Brazil), Nasdaq Composite Indez (NDQ, US), S&P S&P/TSX Composite Inder: (TSX, Canada). In particular, we analyze the scries of markets risk, as mcasured by conditional volatility. Hence, the time series of conditional volatility represent our input data set Y. Information concerning data source and pre-processing given in Appendix H
16

---

Published in Transactions on Machine Learning Research (10/2023)
Regularisation /Fitting Loss
Regularisation /Fitzing Loss
yel
30
yel
20-
20
10-
10-
10
10
A
(a) SS-CASTLE
(b) MS-CASTLE
Figure 8: Behaviour of regularization to fitting loss ratio along À for SS-CASTLE (left) and MS-CASTLE (right), where the x axis is given in log scale.
Methodological approach. As previously noted, owing to the non-convex nature of Problem ( both SS-CASTLE and MS-CASTLE can converge to stationary points that may potentially differ from each other parameter X. Consequently, in our analysis, we seek solutions
across various settings of the for SSCG and MSCG that cxhibit a high degree of persistence across diverse À values in order to mitigate this uncertainty. To achieve this goal, we employ the methodological approach introduced in Section 4.3.2 which proves particularly valuable when no additional data. sets are available to guide the selection of the value of A. Figure Showw shows the behavior of the regularization to fitting loss ratio with respect to X, considering both SS-CASTLE (left) and MS-CASTLE (right). As a meaningful range, we select the values that return a ratio from 0.1 to 1. In this way, we track the change in causal connections when the sparsity-induccis regularization term becomes as important as the model-fitting term. This choice is also due to the fact that conditional volatility series are non-Gaussian and that N As T, similarly to the setting (N, = (10, 128) analyzed in Section 4.3.2) Then, from Figure 8 we select (i) À e [0.004, 0.04] for SS-CASTLE; and (ii) A E [0.003, 0.03] for MS-CASTLE. For cach interval, we pick 10 values for X, and we retrieve the networks of highly persistent edges having Pij 0.95, which show stable sign of the corresponding causal coefficient for all values of A. These edges constitute the causal structures illustrated in the sequel.
5.1 Preliminary Discussion
Before proceeding, let us consider again the causal graphs depicted in Figure I Analyzing the data set Y with SS-CASTLE entails, first and foremost, representing the values of conditional volatility series on a network without edges. Each node of the graph is associated with the value of the i-th risk series, where E [N], taken on specific day i E [T]. The application of SS-CASTLE will yield the SSCG, in which causal edges will convey interactions occurring on a daily basis, possibly with different time lags. However, by examining the system solely in the time domain, won't be able to discern whether causal interactions stem from short-, medium-; or long-term dynamics.
Consider, for example, balance sheet data of listed companies that belong to the stock indices taken into account in our analysis. These data are considered in several asset pricing models, such as those proposed by Fama French (1993); Carhart (1977); Fama & French| (2015), which describe the equity market using lincar models. On one hand, events like carnings relcases, news about company performance, and other firm-specific factors can lead to rapid changes in the short term, causing increased volatility in the related index. On the other hand, fundamental aspects of company, such earnings growth and financial ratios (which are also influenced by macrocconomic conditions), can influence stock prices over longer periods, contributing to valuation trends and implying slower dynamics that evolve at lower frequencies.
In contrast to SS- CASCA,,LE, MS- CASTLE naturally allows us to discriminate between short-, medium-, and long-term interactions through its multiscale analysis using wavelet transform. In this case, the datasct Y will be represented on a network without edges, where the number of nodes is D N TT, with D representing the number of temporal resolutions considered. For example, with reference to the finest time scale d = 1,
17

---

Published in Transactions on Machine Learning Research (10/2023)
the value of node on this graph page will be associated with variations between two consecutive days f and & in the i-th conditional volatility series. This will provide insights into short-term dynamics occurring on a daily temporal resolution. Conversely, on the graph page associated with the time scale d = 3, the value of a node will be linked to medium-term variations between two successive weeks (precisely 4 business days) in the i-th conditional volatility series. Consequently, the causal relationships within the MSCG learned by MS-CASTLE represent causal connections at different time resolutions and time lags.
The analysis presented in the sequel pertains to very specific time period marked by the outbreak and subsequent spread of covid-19, in which it seems natural to expect multiscale causal relationships. Thus, the application of MS-CASTLE seems more appropriate to gain richer insights concerning the dynamics under- lying the system at hand. In particular, short-scale causal interactions will be associated with dependencies between short-term volatility of indices, possibly due to news regarding the contagion or the introduction of new restrictions aimed at containing the spread of the virus. In contrast, medium to long-term causal rela- tionships will intercept dependencies between the riskiness of equity markets tied to pandemic effects, such worsening macroeconomic conditions the impact of mobility restrictions on industrial sectors related to energy (e.g., natural gas and crude oil), material extraction, transportation, and vehicle production.
5.2 Single-scale Causal Analysis
Single
Scale
Scale
511
L
T5X
SIC
SEC
DKC.
TSX.
T5A
(a)
(b)
Figure 9: Highly persistent causal matrices for the single-scale and multiscale where 0.05. Red entries represent positive coefficients, whereas blue ones are negative.
To compare single and multiscale approaches, estimate the causal matrices in both Equations and (4] by using SS-CASTLE and MS-CASTLE respectively. Here, we focus on learning causal graphs from time series using the aforementioned SS-CASTLE method. Figure shows the signed causal matrix made up of persistent coefficients, where for readability reasons we only report the case 0.05. Additional results are provided Appendix T.
In particular, in the matrix shown in Figure 9a the rows represent the parents (sorted according the timestamp), whereas the columns refer to the caused nodes. In this case, based the BIC criterion, set = The upper block of the matrix in Figure 9a concerns lagged causal interactions, while the lower is related to instantancous causal effects. In our study we observe that ares associated with negative causal relations get a weight lower than 0.1 (in module); whereas, most of the surviving edges correspond to autore- gressive causal effects (see also Appendix I). In addition, we find denser causal connections among European and Asia-Pacific countries. Overall, from the single-scale analysis, we cannot find nodes representing major risk drivers within the network, i.c., the considered equity markets show a similar number of outgoing ares (i.e., out-degree).
18

---

Published in Transactions on Machine Learning Research (10/2023)
5.3 Multiscale Causal Analysis
We now focus on multiscale causal analysis, where we analyze the first four temporal resolutions, in accor- dance with the length of our data set, i.c., T 336 obscrvations (see Section 3). As shown by the results in Figure 4 the performance of MS-CASTLE turns out to be robust to different choices of the wavelet. Nevertheless,
consciously select the filter. In general, as shown by Gençay et al. (2010), the choice of the wavelet depends on the length of data, the complexity of the spectral density function, and most impor- tantly the underlying features of data. Hence, similarly to Ren et al. (2021) and references therein, we use Daubechies least asymmetric wavelet with a filter length equal to 8 (db4). Indeed, the reasonable length of this wavelet enhances smoothness (reducing boundary effects in the wavelet coefficients) and allows to deal with even more dynamic spectral density functions. In addition, in this context (and more generally when dealing with financial time series), asymmetric filters are useful to capture stylized features of data and deal with skewed distributions. Here, for instance, the usage of an asymmetric filter is useful to capture rapid increases in market volatility that are followed by gradual declines.
Figure 9b illustrates the highly persistent multiscale causal matrices obtained using the proposed MS- CASTLE method, for the same value of č analyzed above. Additional results are provided in Appendix For the sake of readability, we show separately the diagonal blocks of W corresponding to different scales, i.e. the only elements of W that can be different from zero in Equation
(since no interaction among scales actually takes place). From Figure 9b. we first notice that causal representations different scales show a diverse level of sparsity and, furthermore, persistent causal relations assume only positive values. More in detail, causal interactions are denser at mid-term scales (i.e., 3 and 4, which correspond to 8-16 and 16-32 days, respectively). On the contrary, causal effects turn out to be not persistent at scale 1, which represents a time resolution of 2-4 days. In summary, we found that: (i) the strongest persistent connections appear at scales 3 and 4; (ii) the majority is lagged; and (iii) the most instantaneous relations are associated with weights lower than 0.05 (in absolute value) (see also (ppendixx Th. In addition, we noticed the following behaviors: i) apart from Australia, Asia-Pacific countries are isolated for 5> 0.05; ii) the markets that drive the risk within the network are Brazil, Canada, and Italy. The latter finding can be understood by looking at the number of nonzero entries per market across columns, representing the out-degree of each node. More in detail, the impacts of Brazil and Canada spread across all geographical arcas, while Italy mainly drives the risk within the Eurozone. Finally, from Figure 9b we can notice how the US displays persistent lagged connections.
5.4 Comparison Between Temporal and Multiscale Analysis
The results presented in Sections 5.2 and 5.3 illustrate that, in the case of complex systems such as financial markets, the temporal and multiscale analysis might lead to very diffcrent conclusions. First of all, the inferred SSCG indicates a persistent causal structure al daily frequency, where the strongest connections are autoregressive lagged causal relations. On the contrary, empowered by information concerning the variation of the original signal at different scales, the MSCG shows that causal structures persist at mid-term scales (i.e., 3 and 4), while at short-term scale causal connections are absent. Thus, can conclude that, in our case, the mere application of multiscale-agnostr model leads to a noisy estimate of the causal structure, in which many of the relationships do not persist when decomposing the signal into different temporal resolutions. Drawing from the preceding discussion in Section 5.1 this suggests that within the studied time frame, causal relationships among global equity markets are intertwined with the effect of the pandemic and the prevailing macrocconomic uncertainty. Also, in MSCG we do not obscrve negative causal coefficients as in the case of SSCG, which are somehow difficult to justify during the considered period, since they indicate that an increase (decrease) in the volatility of certain equity market causes a decrease (increase) in that of another market.
Finally, and most importantly multiscale causal analysis allows us to identify the major risk drivers within the network of equity markets during covid-19 pandemic, i.e., Brazil, Canada, and Italy. Interestingly, the
US stock market shows only an impact on Australia, together with an autoregressive effect. In particular, the importance of Canada within the network of stock markets has been underlined by (Ren et al. 2021) well, who conducted study in terms of partial correlation networks. However, since we deal with
19

---

Published in Transactions on Machine Learning Research (10/2023)
causation, result has stronger implications with respect to the aforementioned work. With regards to Brazil, we see that the corresponding stock index shows the highest volatility (see Appendix H), and that its strongest connections (greater than 0.1) are within the American area (see Appendix 1j. We interpret these findings in the context of Canada and Brazil being significant commodity exporters. The pandemic- induced shifts in commodity prices have reverberated through their economies and, consequently, their equity markets. These repercussions can extend to other equity markets, particularly those with robust trade connections. For instance, Canada and the United States have historically demonstrated extensive economic interdependencies, underscored by numerous Canadian firms being cross-listed on U.S. stock exchanges. In a similar vein, Brazil's role as a prominent exporter of agricultural products and machinery to the United States further underscores the potential impact on cross-bordet equity dynamics. Last, but not least, Italy has a high impact within the European area. This nation was among the hardest hit by the pandemic, particularly during its early stages, within Europe. Its key trade relationships encompass a range of sectors, including food, machinery, and vehicles, with major partners like Germany and Francc.
6 Conclusions and Future Research Directions
In this paper, we have proposed novel method to estimate the structure of linear causal relationships at different time scales. By relying upon wavelet transform and non-convex optimization, MS-CASTLE takes explicitly into consideration behaviors of the system at hand spanning diverse time resolutions. Differently from existing causal inference methods, MS-CASTLE looks for linear causal relationships among variations of input signals within multiple frequency bands, and across different time lags. We illustrate that the multiscale-agnoste version of MS-CASTLE, named SS-CASTLE, improves in terms of computational effi- ciency, performance, and robustness over the state of the art. In addition, experiments on synthetic data show that the performance of MS-CASTLE increases with sample availability and is robust to the network size, the noise distribution, and the choice of the wavelets.
The study of the risk of 15 global equity markets, during covid-19 pandemic, shows that USSCASTLE is able to provide useful information about the scales at which causal interactions occur (mid-term scales) and to identify major risk drivers within the system (Brazil, Canada, and Italy). We highlight that the obtained results must be framed in the period of the coronavirus outbreak. Our choice was conscious: given the nonstationary nature of financial markets (Schmitt et al.) 2013), we focused on a narrow period dominated by the pandemic emergency. Thus, the use of different time windows may lead to the estimation of different multiscale causal structure. This observation highlights the need to work on the devclopment of causal inference algorithus capable of handling both the iscale nature of the analyzed system and the nonstationarity of the underlying causal structure (D'A Acunto et al., 2021). In this context, the application of Gaussian processes to model the time dependence of the causal structure has led to some advances (Huang et al. 2015
In addition, the proposed model considers only linear causal relationships generalization to nonlinear inter- actions represents further future work. Here, kernel methods (Shen et al. 2016) and more recently non-linear ICA Monti et al. 2020) have been used to tackle the estimation task. However, previous works only refer to the single-scale case. Also, in this work, we did not consider possible inter-scale cause-effect mechanisms. However, we do not exclude that behaviors of signals at higher frequencies may impact those at lower fre- quencies and vice versa. So, investigating the existence of such causal relationships represents an interesting future research direction. Additionally, would be useful to improve our method for managing multiscale non-stationary causal dynamics. Indcod, it common in different application domains, such as finance, neuroscience, and climatology, to deal with data featured by non-stationarity.
Regards the identifiability, since in the single-scale case Equation {4] coincides with Equation the results achieved for the SVARM hold. Specifically, the causal structure is identifiable in the case of non- Can-Ganssian noise (Shimizu et al. 2011), Gaussian noise with equal or known variances (Peters & Bühhmaann] 201 4 Loh & Bühhlmann 2014), and Gaussian noise with unknown, homogeneous, and heterogeneous variances (Pakk & Kim 2020). In the multiscale case, due to the independence of the frequency bands, we further assume that the underlying contributions coming from different time scales in Equation (4) are observable, the multiscale causal graph is identifiable in all previous cases. However, this is strong assumption that does not hold
20

---

Published in Transactions on Machine Learning Research (10/2023)
in reality. Furthermore, the possible presence of serial correlation in the noise decomposition hinders the identifiability of a multiscale causal model. Hence, since we did not establish any identifiability result in the case of unobservable contributions and multiscale noise with serial correlation, we plan to investigate this intriguing point in future work.
Related to the identifiability issue, the guarantee recovering the underlying causal structure provided by least-squares- based methods relates assumptions on the noise variances (Theorems 7 and 9 in Loh & Bühhmannn 2014h. As recently highlighted by Ng et al. (2023), data standardization facilitates the violation of such assumptions. Hence, in the case of standardized data the proposed methods (and more generally least-squares-basee methods) do not have any guarantec of recovering the true causal structure.
Finally, the results of the case study show how MS-CASTLE can be used to support portfolio risk manage- ment. Indeed, depending on their investment horizon, investors could use the proposed methodology to make risk-aware decisions regarding their portfolios, from causal perspective and without any prior assumption about the scale of analysis.
References
Mark Alber, Adrian Buganza Tepole, William R Cannon, Suvranu De, Salvador Dura-Bernal, Krishna Garikipati, George Karniadakis, William W Lytton, Paris Perdikaris, Linda Petzold, et al. Integrating machine learning and multiscale modeling perspectives, challenges, and opportunities in the biological, biomedical, and behavioral sciences. NPJ Digital Medicine, 2(1):1-11, 2019.
Michel Besserve, Bernhard Schölkopf, Nikos K Logothetis, and Stefano Panzeri. Causal relationships between frequency bands of extracellular signals in visual cortex revealed by an information theoretic analysis. Journal of Computational Neuroscience, 29:547-566, 2010.
Monica Billio, Mila Getmansky, Andrew W Lo, and Loriana. Pelizzon. Econometric measures of connectedness and systemic risk in the finance and insurance sectors Journal of Financial Economics, 104(3):533--55) 2012.
Tim Bollerslev. Generalized autoregressive conditional heteroskedasticity. Journal of Econometrics, 31 31)): 307 327, 1986.
Stephen Boyd, Neal Parikh, and Eric Chu. Distributed Optimization and Statistical Learning via the Alter- nating Direction Method of Multipliers. Now Publishers Inc, 2011.
Peter Bühhlmann, Jonas Peters, Jan Ernest, et al. CAM: Causal additive models, high-dimensional order search and penalized regression. The Annols of Statistics, 42(6):2526 2556, 2014.
Richard Byrd, Peihuang Lu, Jorge Nocedal, and Ciyou Zhu. A limited memory algorithm for bound constrained optimization. SIAM Journal on Scientific Computing, 16(5):1190-120) 1995.
Laurent Calvet and Adlai Fisher. Forccasting multifractal volatility. Journal of Econometrics, 105(1):27-58, 2001.
Mark M Carhart. On persistence in mutual fund performance. The Journal of Finance, 52(1))57 82, 1997.
David Maxwell Chickering. Optimal structure identification with greedy search. Journal of Machine Learning Research, 35055,5 2002.
Max Chickering, David Heckerman, and Chris Meek. Large-sample learning of Bayesian networks is NP-hard. Journal of Machine Learning Research, 5, 2004.
Philippe Ciuciu, Patrice Abry, and Biyu J He. Interplay between functional connectivity and scale-free dynamics in intrinsic IMRI networks. Neuroimage, 95:248 263, 2014.
Gabricle D'Acunto, Paolo Bajardi, Francesco Bonchi, and Gianmarco De Francisci Morales. The evolving causal structure of equity risk factors. In Proceedings of the Second ACM International Conference on Al in Finance, pp. 1-8, 2021.
21

---

Published in Transactions on Machine Learning Research (10/2023)
Gabriele D'Acunto, Gianmarco De Francisci Morales, Paolo Bajardi, and Francesco Bonchi. Multiscale non-stationary causal structure learning from time series data. ar arXiv preprint arXin:2208.. 14989, 2022.
Eugene F Fama and Kenneth R. French. Common risk factors in the returns on stocks and bonds. Journal of Financial Economics, 33(1):3-56 1993.
Eugene F Fama and Kenneth R French. A five-factor asset pricing model. Journal of Financial Economics, 116(1):1 22, 2015.
Ramazan Gençay, Nikola Gradojevic, Faruk Sclçuk, and Brandon Whitcher. Asymmetry of information flow between volatilities across time scales. Quantitative Finance, 10(8):895-99) 915, 2010.
John Geweke. Measurement of linear dependence and leedback between multiple time series. Journal of the American Statistical Association, 77(378):304 -33,3) 1982.
John F Geweke. Measures of conditional linear dependence and feedback between time series. Journal of the American Statistical Association, 993399))0 1984.
Clark Glymour, Kun Zhang, and Peter Spirtes. Review of causal discovery methods based on graphical models. Frontiers in Genetics, 10:524, 2019.
Donald Goldfarb, Shiqian Ma, and Katya Scheinberg. Fast alternating linearization methods for minimizing the sum of two convex functions. Mathematical Programming, 141(1):349 (1):4444382)) 382, 2013.
Mingming Gong, Kun Zhang, Bernhard Schölkopf, Dacheng Tao, and Philipp Geiger. Discovering temporal causal relations from subsamplod data. In International Conference Machine Learning, pp. 1898-1906.
PMLR, 2015.
Clive WJ Granger. Investigating causal relations by econometric models and cross-spectral mots Econo- metrica: Journal of the Econometric Society, pp. 424 438, 1969.
Shuixia Guo, Anil K Seth, Keith M Kendrick, Cong Zhou, and Jianfeng Feng. Partial Granger causal- ity Ellminata exogenous inputs and latent variables. Journal of Neuroscience Methods, 172(1):79 93, 2008.
Kazuharu Harada and Hironori Fujisawa. Sparse estimation linear non-Gaussian acyclic model for causal discovery. Neurocomputing, 459:223 233, 2021.
David Heckerman, Dan Geiger, and David M Chickering. Learning Bayesian networks: The combination of knowledge and statistical data. Machine Learning, 20(3):197-2 14319977.,,,.. 1995.
Johan Himberg, Aapo Hyvärinen, and Fabrizio Esposito. Validating the independent components of neu- roimaging time series via clustering and visualization. Neuroimage, 2213) 2004.
Patrik Hoyer, Dominik Janzing, Joris M Mooij, Jonas Peters, and Bernhard Schölkopf. Nonlinear causal discovery with additive noise models. Advances in Neural Information Processing Systems, 21:689-696, 2008.
Biwei Huang, Kun Zhang, and Bernhard Schölkopf. Identification of time-dependent causal model: A Gaussian process treatment. In Tuvetty-Fourth International Joint Conference on Artificial Intelligence, 2015.
Biwei Huang, Kun Zhang, Yizhu Lin, Bernhard Schölkopf, and Clark Glymour. Generalized score functions for causal discovery. In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pp. 1551 1560, 2018
Biwei Huang, Kun Zhang, Jiji Zhang, Joseph Ramsoy, Ruben Sanchez-Romero, Clark Glymour, and Bern- hard Schölkopf. Causal discovery from naterogenoons/noossatiiono data. Journal of Machine Learning Research, 21 (89):1 109::1-5,,, 2020.
22

---

Published in Transactions on Machine Learning Research (10/2023)
Aapo Hyvarinen. Fast and robust fixed-point algorithms for independent component analysis. IEEE Trans- actions on Neural Networks, 10(3):626 634, 1999.
Aapo Hyvärinen, Kun Zhang, Shohei Shimizu, and Patrik O Hoyer. Estimation of a structural vector autoregression model using non-Gaussianity. Journal of Machine Learning Reseorch, (1)),, 2010.
Jaime Ide, Fábio Cappabianco, Fabio Faria, and Chiang-shan R Li. Detrended partial cross correlation for brain connectivity analysis. Advances Neural Information Processing Systems, 30, 2017.
Joshua Jacobs, Michael J Kahana, Arne D Ekstrom, and Itzhak Fried. Brain oscillations control timing of 2007.
single-ncuron activity in humans. Journal of Neuroscience,
Steve Kalke and W-D Richter. Simulation of the p-generalized Gaussian distribution. Journal of Statistical Computation and Simulation, (34))6641-667, 2013.
Marcin Jan Kaminski and Katarzyna J Blinowska. A new method of the description of the information flow in the brain structures. Biological Cybernetics, 65(3):203-2N) 210, 1991.
Rabeh Khalfaoui, Mohamed Boutahar, and Heni Boubaker. Analyzing volatility spillovers and hedging between oil and stock markets: Evidence from wavelet analysis. Energy Economics, 19:540-549,, 2015.
Lutz Kilian and Helmut Lütkepohl. Structural vector autoregressive analysis. Cambridge University Press, 2017.
Jarostaw Kwapień and Stanislaw Droidž. Physical approach to complex systems. Physics Reports, 515(3-4): 115-226, 2012.
Sébastien Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon Lacoste-Julien. Gradient-based neu- ral DAG learning. In International Conference Learning Representations, 2020. URL https:
//openreview. www..e//ooormmdddriibee.ee
Lixia Loh. Co-movement of Asia-Pacific with European and US stock market returns: A cross-time-rrequeecy analysis. Research in International Business and Finance, 29:1-13, 2013.
Po-Ling Loh and Peter Bühlmann. High-dimenssonall learning of linear causal networks via inverse covariance cstimation. Journal of Machine Learning Research, 15()):3065-3105, 2014.
Rosario N Mantegna and H Eugene Stanley. Introduction to econophysies: Correlations and complezity in finance. Cambridge University Press, 1999.
Alessio Moneta, Doris Entner, Patrik O Hoyer, and Alex Coad. Causal inference by independent component analysis: Theory and applications. Oxford Bulletin of Economics and Statistics, 7B(5):705-733, 2013.
Alessio Moneta, Gianluca Pallante, al. Identification of structural VAR models via independent component analysis: A performance evaluation study. Technical report, Laboratory of Economics and Management (LEM), Sant' Anna School of Advanced Studies, 2020.
Ricardo Pio Monti, Kun Zhang, and Aapo Hyvärinen. Causal discovery with general non-linear relationships using non-linear ICA. In Uncertainty in Artificial Intelligence, pp. 186 195. PMLR, 2020.
Guy P Nason. Wavelet methods in statistics with R. Springer, 2008
Guy P Nason and Bernard W Silverman. The stationary wavelet transform and some statistical applications.
In Wavelets and statistics, pp. 281-299. Springer, 1995.
Tgnavier Ng and Kun Zhang. Towards federated Bayesian network structure learning with continuous op- timization. In International Conference on Artificial Intelligence and Statistics, pp. 8095-8111. PMLR, 2022.
Ignavier Ng, AmirEmad Ghassami, and Kun Zhang. On the role of sparsity and DAG constraints for learning linear DAGs. Advances in Neural Information Processing Systems, 33:17943- 17954, 2020.
23

---

Published in Transactions on Machine Learning Research (10/2023)
Ignavier Ng. Biwei Huang, and Kun Zhang. Structure learning with continuous optimization: A sober look and beyond. XrXie preprint or 230:2334002144 2023.
Guido Nolte, Andreas Ziehe, Vadim V Nikulin, Alois Schlögl, Nicole Krämer, Tom Brismar, and Gelaus-Robers Müller. Robustly estimating the flow direction of information complex physical systems. Physical Review Letters, 103423): 2008.
Roxana. Pamfil, Nisara Sriwattanaworschai, Shaan Desai, Philip Pilgerstorfer, Konstantinos Georgatzis, Paul Beaumont, and Bryon Aragam. DYNOTEARS: Structure learning from time-series data. In International Conference on Artificial Intelligence and Statistics, pp. 1595-1605. PMLR, 2020.
Gunwoong Park and Youngwhan Kim. Identifiability of Gaussian linear structural equation models with homogeneous and heterogeneous error variances. Journal of the Korean Statistical Society, 49:276 292, 2020.
Judea Pearl. Causality. Cambridge University Press, 2009.
Grace CY Peng, Mark Alber, Adrian Buganza Tepole, William R Cannon, Suvranu De, Savador Dura-Ber Bara-Bernal, Krishna Garikipati, George Karniadakis, William W Lytton, Paris Perdikaris, al. Multiscale modeling meets machine learning: What can we learn? Archives of Computational Methods in Engineering, 28(3):
1017-1037, 2021.
Donald B Percival and Harold O Mofjeld. Analysis of subtidal coastal sea level fluctuations using wavelets. Journal of the American Statistical Association, 92(439):868 880, 1997.
Donald B Percival and Andrew T Walden. Wavelet Methods for Time Series Analysis, volume 4. Cambridge University Press, 2000.
Jonas Peters and Peter Bühlmann. Identifiability of Gaussian structural equation models with equal error variances. Biometrika,
2014.
Jonas Peters, Joris M Mooij, Dominik Janzing, and Bernhard Schölkopf. Causal discovery with continuous additive noise models. Journal of Machine Learning Research, 2014.
Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Elements of causal inference: Foundations and learning algorithms. The MIT Press, 2017.
Yinghua Ren, Wanru Zhao, Wanhai You, and Kaikai Zhai. Multiscale and partial correlation networks anal- ysis of risk connectedness in global equity markets. Physica A: Statistical Mechonics and its Applications, 573:125911, 2021.
Jakob Runge. Discovering contemporaneous and lagged causal relations in autocorrelated nonlinear time series datasets. In Conference on Uncertainty in Artificial Intelligence, pp. 139871397. PMLR, 2020.
Jakob Runge, Sebastian Bathiany, Erik Bollt, Gustau Hamps-Valls, Dim Coumou, Ethan Deyle, Clark Glymour, Marlene Kretschmer, Miguel D Mahocha, Jordi Muñoz-Marí, et al. Inferring causation from time series in Earth system sciences. Nature Communécations, 10(1):1-13, 2019a.
Jakob Runge, Peer Nowack, Marlene Kretschmer, Seth Flaxman, and Dino Sejdinovie. Delecting and quanti- fying causal associations in large nonlinear time series datasets. Science Advances, 5(11)rea 2019b.
Thilo A Schmitt, Desislava Chetalova, Rudi Schäfer, and Thomas Guhr. Non-stationarity in financial time scries: Generic features and tail behavior. EPL (Europhysics Letters), 103(5):58003, 2013.
Bernhard Schölkopf, Francesco Locatello, Stefan Bauer, Nan Rosemary Ke, Nal Kalchbrenner, Anirudh Goyal, and Yoshua Bengio. Toward causal representation learning. Proceedings of the IEEE, 109(5): 2021.
Gideon Schwarz. Estimating the dimension of model. The Annals of Statistics, pp. 461-464, 1978.
24

---

Published in Transactions on Machine Learning Research (10/2023)
Yanning Shen, Brian Baingana, and Georgios B Giannakis. Nonlinear structural vector autoregressive models for inferring effective brain network connectivity. ar preprint 1et166600655,, 2016.
Shohei Shimizu, Patrik O. Hoyer, Aapo Hyvärinen, and Antti Kerminen. A linear non- non-Gausiian acyclic December 2006. ISSN
model for causal discovery. Journal of Machine Learning Research,
1532-4435.
Shohei Shimizu, Takanori Inazumi, Yasuhiro Sogawa, Aapo Hyvärinen, Yoshinobu Kawahara, Takashi Washio, Patrik O. Hoyer, and Kenneth Bollen. DirectLINGAM: A direct method for learning a lin- non-Gaussian structural equation model. Journal of Machine Learning Research, 12:1225- 141122--248,, 1248, July
2011. ISSN 1532-4435.
Christopher A Sims. Macroeconomics and reality. Econometrica: Journal of the Econometric Society, pp. 1-48, 1980.
Peter Spirtes, Clark N Glymour, Richard Scheines, and David Heckerman. Causation, Prediction, and Search. MIT Press, 2000.
Cornelis J Stam and Elisabeth CW van Straaten. Go with the flow: Use of directed phase lag index (dPLI) to characterize patterns of phase relations in large-scale model of brain dynamics. Neuroimage, 62(3): 1415-1428, 2012.
Matthew Vowels, Necati Cihan Camgoz, and Richard Bowden. D'ya like DAGs? A survey on structure learning and causal discovery. ACM Computing Surveys, 55(4):1-36, 2022.
Gang-Jin Wang, Chi Xie, Min Lin, and H Eugene Stanley, Stock market contagion during the global financial crisis: multiscale approach. Finance Research Letters, 22:163 168, 2017.
Junfeng Yang and Xiaoming Yuan. Lincarized augmented Lagrangian and alternating direction methods for nuclear norm minimization. Mathematics of Computation, 82(281):301 329, 2013.
Xing Yang, Chen Zhang, and Baihua Zheng. Segment-wise time-varying dynamic Bayesian network with graph regularization. ACM Transactions Knowledge Discovery from Data (TKDD), 16(6):1-23, 2022.
Yue Yu, Jie Chen, Tian Gao, and Mo Yu. DAG-GNN: DAG structure learning with graph neural networks. In International Conference on Machine Learning, pp. 7154-7163. PMLR, 2019.
Xun Zheng, Bryon Aragam, Pradeep Ravikumar, and Eric P Xing, DAGs with no tears: Continuous optimization for structure learning. In Advances in Neural Information Processing Systems, volume 31, 2018.
Xun Zheng, Chen Dan, Bryon Aragam, Pradeep Ravikumar, and Eric Xing. Learning sparse nonparametric DAGs. In International Conference on Artificial Intelligence and Statistics, pp. 3414-3425. PMLR, 2020.
Shengyu Zhu, Ignavier Ng, and Zhitang Chen. Causal discovery with reinforcement learning. In International Conference Learning Representations, 2020.
25

---

Published in Transactions on Machine Learning Research (10/2023)
Appendix A The Wavelet Transform
The discrete wavelet transform (DWT, Percival & Mofjeld 1997 is mathematical technique that projects a time series onto a set of orthonormal basis functions, known as wavelets. This projection results in a series of wavelet coefficients that capture the information from the time series across different frequencies and time intervals.
Unlike the discrete Fourier transform (DFT), which uses sine and cosine functions for data projection and yields coefficients that not time-varying, the DWT benefits from time resolution by employing basis functions having compact support, i.e, function oscillating over brief time interval.
As an example, consider the Haar wavelet, having filter h = (5 - -/22 of length K = 2. It is straightforward to see that this filter (i) has zero-sum Σx h4 = O; (ii) has unit energy hh" = 1, and (iii) is orthogonal to its even shifts Σx h1/1:220 n E Z, n + O.
These properties are very important since ensure that (i) the wavelet coefficients identify changes in the data, (ii) the transformation preserves the input signal variance without neither adding nor suppressing information, and (iii) the DWT enables the analysis of the signal at multiple time resolution by decomposing it on orthonormal basis functions.
Associated with the decomposition filter, there is the scaling filter g = ( ) that similarly to h preserves data variance and is orthonormal. However, if h can be thought of as a local differentiating operator, g can be interpreted as a local averaging operator.
The DWT applies and g in a series of filtering and downsampling operations to the input time series having a dyadic length T = 2D, with E N. These operations yield D wavelet coefficients, effectively decomposing the original time scries into components living in both time and frequency domains.
The two filters h and g are also known as high-pass and low-pass filters. This is due to the fact that their transfer functions, i.e,, their DFT, favor high and low frequencies, respectively. For instance, h keeps high frequencies and suppresses low frequencies. Hence, by being applied together, h and g decompose the input time series in high- and low-frequency representations, preserving all its content.
In order to understand the filtering and downsampling operations, let us consider the sequence s = (1, 4) of length T 4 and the aforementioned Haar wavelet. As a first step, the filtering operation Specifically, starting from 8, we obtain the vector of wavelet coefficients at scale d = 1 by convolving (with non-overlapping windows) h with s, b (b), bb) = (- 3 v/), where bj
applied.
1 1 1 2 and v2
v2
v2
b/= = 1 3- -2 4. At this point, we apply the downsampling operation by convolving in the same way g with s, obtaining the first vector of scaling coefficients, c² = (cl, c)) = 3
v2
( 5 We proceed in the same way for the subsequent scale level d = 2. These two operations can be compactly described by the following formulas:
K 1/4
K-1
hach 22 mod
4=
Σ
gh 26-2111 mod
d-1
k= 0
k=
where c and E [ -1o
1)o-
It is easy to see from the previous example that DWT is sensitive to the origin, i.e., shift in would lead to different wavelet coefficients. Morcover, a question arises as to whether exclusively focusing on even differences might overlook certain patterns within the data. Starting from the DWT, the stationary wavelet transform (SWT, Nason & Silverman 1995) builds on the basic observation that, in some cases, it is useful to retain both odd and even differences at each scale. In addition, contrast to DWT, SWT is not origin- sensitive, i.e., translation invariant, and avoids decimating the input signal at each scale, i.c., halving its length at each iteration. We refer the interested reader to Chapter 2.9 in Nason (2008) for further details.
26

---

Published in Transactions on Machine Learning Research (10/2023)
Appendix B Comparison with DYNOTEARS
Figure depictss the swarm plot concerning the number of edges and the SHD (structural Hamming distance)) of the estimated matrices of causal coefficients. 'he latter metric indicates the number of modifications needed to retrieve the ground truth from the estimated causal graph (the lower, the better). First, we observe that both models converge to causal networks of similar size. In addition, by looking at SHD, we notice that dagness function linearization ion does not cause worsening in estimation accuracy.
01 00tee/00
10
100
10
30
100
Nodes
Nodes
(a) W
10
90
100
10
30
100
Nodes
Nodes
(b) W
Figure 10: Swarm plots regarding the number of edges (left) and SHD (right) of the estimated causal matrices.
411
-31
Misinn
Reverse
Extra
101
10
01
20
90
30
100
Nodes
Nodes
Nodes
(a) W
PINOTEANS
Micsing
Extra
Nodes
(b) W
Figure 11: Swarm plots regarding the building blocks of SHD (from the left: extra, missing and reverse edges wrt the ground truth) associated with the estimated causal matrices wo (a) and W1 (b). Please notice that in the case of lagged causal interactions, we cannot observe reverse edges.
The comparison in terms of accuracy between DYNOTEARS and SS-CASTLE, is further detailed below. In particular, Figure II shows the contribution of eztra, missing, and reverse edges the SHD. More precisely:
27

---

Published in Transactions on Machine Learning Research (10/2023)
eatra edges are estimated edges not encompassed in the causal graph skeleton;
missing edges are causal connections present in the ground truth that have not been retrieved, not even with a wrong direction;
reverse edges are those connections estimated with a wrong direction.
Please, refer to Appendix C for detailed definition of the considered metrics. We observe that the number of missing edges by far the most dominant component (especially in the case of larger networks). Overall, the models perform similarly across the three components. It is worth noticing that, since W is acyelic by definition, we cannot have reverse edges.
Appendix C Metrics
We evaluate the learned causal graphs in terms of F1-score and structural Hamming distance (SHD), com- puted from the adjaconcy matrices of the learned causal graphs. Let us consider
the condition positives (p) as the number of edges present in the ground truth;
the condition negatives (n) as the number of ares that are not in the ground truth;
the number of estimated edges (nnz);
the true positives (tp) the number of estimated edges present in the ground truth and having the
correct direction;
the reversed (r) as the number of learned edges present in the ground truth but with the opposite direction;
the false positives (fp) as the number of learned extra edges that are not in the undirected skeleton of the ground truth;
the missing edges (c) as the number of edges of the skeleton of the learned graph that extra from
the skeleton of the ground truth;
the extra edges (m) as the number of edges of the skeleton of the learned graph that are missing from the skeleton of the ground truth.
At this point, define the false discovery rate as fdr = +r++ inn, the true positive rate as tpr tp/p, and the false positive rate for = fp)nn. Hence, the F1-score is equal to 2 (1-fdr)tpr and the SHD is equal to f-fdr+tpr
+ m
In practical terms, to compute the previous metrics we proceed as follows. Starting from the estimated W we remove the zero off dif-diagnaal blocks that come into play when we build the diagonla matrices W (see Section 2.1 In this way, we consider only those wi blocks that can contain entries different from zero, corresponding to ares in the causal graph. At the end of this step, end up with reshaped version of W, having size Px (X+1) Then, we transform the latter 4-dimensional tensor in logic one, where entry is equal to 1 iff it corresponds to a nonzero entry in the original tensor, and 0 otherwise. In this way, we retain only the structural information, that can be interpreted as the multiscale counterpart of the adjacency. We repeat the same steps for the ground truth.
Then, is possible count the quantities listed above, needed for the computation of the Fl-score and SHD, by cycling over the time scales and time lags.
28

---

Published in Transactions on Machine Learning Research (10/2023)
Appendix D
Hyper-parameters
Table provides the chosen values for À during the experimental assessment of SS-CASTLE. In addition, Table shows the values of À used in the experimental assessment of MS-CASTLE. In our experiments, we considered E (0.0001, 0.0005, 0.001, 0.005, 0.01, 0.05, 0.1.
A
Nodes
T- 100
T 500
T-1000
1.0
0.50
0.10
0.10
1.5
0.10
0.05
0.05
10
2.0
0.50
0.10
0.10
2.5
0.10
0.05
0.01
100.0
0.10
0.05
0.01
1.0
0.50
0.10
0.10
1.5
0.10
0.05
0.05
2.0
0.50
0.10
0.10
2.5
0.10
0.05
0.01
100.0
0.10
0.05
0.01
1.0
0.50
0.10
0.10
1.5
0.50
0.05
0.05
50
2.0
0.50
0.10
0.10
2.5
0.10
0.05
0.05
100.0
0.10
0.05
0.01
Table Selected values for X for each of the considered parameters combinations (7, N,,N,D p).-
29

---

Published in Transactions on Machine Learning Research (10/2023)
A
Nodes
p
Transform
T-128
T-512
T-8192
None
0.0100
0.0050
0.0050
swt-db3
0.0500
0.0100
0.0100
1.0
swt-db5
0.0500
0.0100
0.0100
swt-db7
0.0500
0.0100
0.0100
5
swt-symā
0.0500
0.0050
0.0100
swt-coif2
0.0500
0.0050
0.0100
None 0.0050
0.0050
0.0010
swt-db3
0.0050
0.0010
0.0010
2.0
swt-db5
0.0050
0.0010
0.0001
swt-db7
0.0050
0.0010
0.0005
swt-sym5
0.0050
0.0010
0.0001
swt-coif2
0.0050
0.0010
0.001
None 0.0500
0.0100
0.0050
swt-db3
0.0500
0.0100
0.0100
1.0
swt-db5 0.0500
0.0500
0.0100
swt-db7 0.0500
0.0500
0.0100
10
swt-sym5
0.0500
0.0500
0.0100
swt-coif2
0.0500
0.0100
0.0100
None 0.0050
0.0010
0.0010
swt-db3
0.0100
0.0050
0.0010
2.0
swt-db5
0.0100
0.0050
0.0010
swt-db7
0.0100
0.0050
0.0010
swl-sym55
0.0100
0.0050
0.0005
swt-coif2
0.0100
0.0050
0.0010
None
0.1000
0.0100
0.010
swt-db3
0.0500
0.0500
0.0100
1.0
swt-db5
0.0500
0.0500
0.0100
swt-db7 0.0500
0.0500
0.0500
30
swt-symö
0.0500
0.0500
0.0100
swl-coif2
0.0500
0.0500
0.0100
None 0.0500
0.0050
0.0005
swt-db3
0.0100
0.0050
0.0005
2.0
swt-db5
0.0100
0.0050
0.0010
swt-db7
0.0100
0.0050
0.0010
swt-sym5
0.0100
0.0050
0.0010
swt-coif2
0.0100
0.0050
0.0005
Table 2: Selected values for for the tested versions MS-CASTLE, as indicated by the column Transform.
30

---

Published in Transactions on Machine Learning Research (10/2023)
Appendix E
Additional Comparison with
Methods
Sensitivity to data set size.
Figure 12 provides additional details concerning the structural mistakes made by the models. In particular, we provide the comparison with VAR-DIrectLINGAM, VAR-ICALINGAM, and PCMCI+ in terms of SHD, and its building blocks, i.e., extra, missing and reverse edges. Overall, SS-CASTLE outperforms the baselines also in terms of the considered structural metrics. Looking at the instantancous causal interactions (top rows), we notice how non-Gaussian methods lend to estimate greater number of extra and reverse edges, even when T' = 1000. Conversely, PCMCI+ exhibits a low number of extra and reverse edges, making it competitive with the proposed method. With regards to missing edges, all the models but PCMCI+ tend to perform similarly as T grows. Indeed, we observe that the main driver in PCMCI+ estimation error is the number of missing edges, which remain high across all the considered settings. Furthermore, the results show that non-Gaussian methods display a dependence on the value of p. In order to better interpret the values of SHD, consider that in case N 30 and .85, average only 65 entries of W are different from zero due to the acyclicity requirement.
Concerning lagged causal interactions (bottom rows), we see that the non-Gaussian methods are prone to return solutions characterized by a large number of extra edges. Furthermore, as is the case with instanta- neous relations, the primary factor determining the SHD shown by PCMCI+ is the number of missing edges. Regarding SS-CASTLE, the number of missing edges turns out be the major contributor to SHD in case of small data sets (T 100). With the increase of T, extra and missing edges start to contribute similarly to the aforementioned structural metric.
Sensitivity to network size.
Figure 13 provides further information regarding the estimated structure when we vary the network size. Also this case we see that SS-CASTLE outperforms the considered basclines. Concerning the learning of the matrix of instantaneous relations, we see that non-Gaussian methods are prone to estimate greater number of extra and reverse edges. As for PCMCI+, the major driver in SHD is the number of missing edges. Specifically, the values shown by the latter method are the worst. With regards SS-CASTLE, the main component of SHD is the number of missing edges. As far as the learning of the matrix of lagged causal interactions, the results show that, even though the non-Gaussian met and SS-CASTLE display a similar number of missing ares, overall SS-CASTLE is more robust to extra edges. Finally, in the case of PCMCI+, the number of missing edges remains the primary source of crror.
Appendix F
Multiscale Causal Structure
The multiscale causal structure determining the data generated for the experimental assessment of MS- CASTLE consists of 4 different time scales. In addition to instantancous causal interactions, to test also the performance of MS-CASTLE in the presence of lagged interactions, we set the autoregressive lag I, = 1, for each time scale.
In detail, over the first time scale W! has entries [Wik.+1 = 0.6 and zero elsewhere, i E [N - 1]; W has entries (i) Wi],, = 0066, with i E [N], (ii) [Wilast1 = 0.3, with j E [N - -1, (iii) [W]] = 0.3, with k e [N 1], and (iv) zero elsewhere.
At the second scale W3 has entries [W 00.3, with E [N-2], and zero elsewhere; Wi has entries (i) [W -0.5, with i [N], (ii) [Wilas+2 = 0.4, with j E [N 2],, (iii) [Wi] -0.4, with A [N 2]. and (iv) zero elsewhere.
Over the third and fourth scale we set W with entries (i) [W]] = 0.5, with E [N], (ii) [Wil = -04, set with j E [N - 1], (ii) [Wi] = 0.4, with & e [N - 1], and (iv) zero elsewhere; Wi with entries (i) [W = 0.7, with E [N], and zero elsewhere.
Finally, in our experiments, we consider N E (5, 10,
31

---

Published in Transactions on Machine Learning Research (10/2023)
7 100
T 500
1.00
100
Valk (heett WGAY
Lak
100
199
100
SHD
1.0
1000
1.5
2.0
25
100.0
10
1.5
103.0
(a) SHD.
7 100
T 500
1.00
SAR: Diee- WGAV ICALIEAAA
100
Exox Lag
100
Exox Lik
106.0
1.5
2.0
255
100.0
10
155
103.0
(b) Extra edges.
T' 100
T 000
2 2000
100-
Masiug Lag
60
100.
Mosing Lag
00-
10
1.5
20
25
1060
10
1.5
2.0
7.5
100
10
2.0
1030
Missing edges.
2003
T 200
T- 1000
30
Merets, Lag
10
WWKDbeetL/MAA VAS ICA, 2MCAM
Roaens, Lag:
POVLL CASTLE
30-
100.0
(d) Reverse edges.
Figure 12: Comparison with the baseline methods in the estimation of the causal matrices n (top) and W
(bottom) in terms of SHD, and its building blocks, i.c., b eextra edges, X missing edges, and i reverse edges.
Each subfigure depicts the case when N 30 and the number of time series varies {100, 500, 1000).
32

---

Published in Transactions on Machine Learning Research (10/2023)
10
N
50
15
VAR CeA
100
210-
Lag
50
300
SHD
200
190
5000
SHD Lag
50
100
2.5
100.0
1.0
100.0
1.0
1.5
2.0
100.0
(a) SHD.
00
100-
WAL-CA MCUMAA
Eem Lz
POWOH
CASTEE
30
20
10
Estre Lag.
200-
1.5
2.0
2.5
100.0
1.0
2.5
100.0
1.0
1.5
2.0
2.5
100.0
(b) Extra edges.
10
N 30
N-
7.5
Mashug Lake
5.0
0.0
Moking Lag
5.0
50
1.55
1.5
2.0
25
100.0
2.5 100.0
103.0
(c) Missing edges.
N
7.5-
freens, lag
0.5-
Moeese Lay
000
0.00
000
1.0
1.5
2/a
9.s
106.0
110
9/5
100.0
2's 1000
(d) Reverse edges.
Figure 13: Comparison with the baseline methods in the estimation of the causal matrices Wo (Lop) and W (bottom) in terms of ofla SHD, and its building blocks, i.e., b eextra edges, I missing edges, and d reverse edges. Each subfigure depicts the case when T 1000 and the number of time scries N varies in (10, 30, 50).
33

---

Published in Transactions on Machine Learning Research (10/2023)
Appendix G MS-CASTLE: Additional Results on Synthetic Data
N
N 10
N
Gairsian care
1.0
0.5
0.0
non Causian c3t
1.0
0.5
MSCASTIE
MECASTLE vl 417
MOLASTLE not 9-1
0.0
MECASTIE sat oxf
128
612 8192
128
512
8102
128
512 8192
T
T
Figure 14: Comparison of different versions of MS-CASTLE in terms of normalized SHD (i.e., SHD divided by the number of edges in the ground truth), obtained in the experimental settings described in Section 4.3.1 Subplots on the top refer to the Gaussian setting, whereas those on the bottom to the non-Gaussian one. Subplots on the left refer to the case N 5, those on the center to N = 10, and those on the right to N 30. Dashed lines represent the inter-quartile range.
34

---

Published in Transactions on Machine Learning Research (10/2023)
2" 225
T=
F 120
00.75
0.75
0.30
(a) Optimal decomposition.
7-
300
7
1.00
EURI
01.00
E000
0..00
(b) coif2.
(c) db3.
0.70
Grasten
8220
8.55
00
117
12
107 10
18
10
10
10
167
(d) db7.
(e) sym5.
Figure 15: Behaviour of F1-score along the considered values of X, obtained in the experimental settings described Section 4.3.1 Subplots on the top refer to the Gaussian setting, whereas those on the bottom the non-Gaus one. Subplots on the left refer to the case T 128, those on the center T 512, and those the right to T 8192. Dashed lines represent the inter-quartile range, while colors refer to different values of N.
35

---

Published in Transactions on Machine Learning Research (10/2023)
Appendix H
Financial Data
140
120
Value
100
80
2020-00
2020-03
2020-05
2020-07
2020-09
2020-11
2021-01
2021-03
2021-05
Ordinaries Index
TAIEX Index
FTSE 100 Index
Hang Seng Index
DAX Index
RTS Index
Nikkei 225 Index
FTSE MIB Index
Bovespa Index
Shanghai Composite Index
IBEX Index
Nasdaq Compusite Index
Straits Times Index
CAC 40 Index
S&P S&P/TSX /SX Composite Index
Figure 16: Behavior of equity markets during the considered period. All indexes are rebased to 100.
We consider data concerning 15 global equity markets at daily frequency. To focus covid-19 pandemic period, we restrict our attention to observations from January 2020, the 2nd to April 2021, the 30th. In our analysis, we deal with the following markets: All Ordinaries Index (AOR, Australia), Hang Seng Index (HSI, Hong Kong), Nikkei 225 Indez (NKX, Japan), Shanghai Composite Inder (SHC, Shanghai), Stroits Times Inder (STI, Singapore), TAIEX Inder (TWSE, Taiwan), DAX Inder (DAX, Germany), FTSE MIB Inder (FMIB, Italy), IBEX Index (IBEX, Spain), CAC 40 Inder (CAC, France), FTSE 100 Inder (UKX, UK), RTS Inder (RTS, Russia), Bovespa Inder (BVP, Brazil), Nasdag Composite Inder (NDQ, US), S&P/TSX Composite Indez (TSX, Canada). The data has been downloaded from Stood Figure 16 depicts the behavior of the indexes during the considered time window. In particular, the indexes plummeted during the first msnths of 2020 and, subsequently, they showed a second downturn during October 2020. In addition, Table 3 provides summary statistics. Overall, according to risk-adjusted return Sortino ratid and average compounded return to max drawdown ratio (ACR/MDD), TWSE and NDQ outperform the rest of the indexes. Moreover, we see that annualized average compounded returns largely vary across the considered instruments: while IBEX and UKX are the worst performing, TWSE and NDQ the most profitable. Furthermore, all indexes show high level of volatility. Among the others, BVP and RTS are the most volatile indexes. Last but not least, all indexes suffer heavy losses during the analysed period, as shown by max drawdown metric (MDD). Interestingly, SHC shows the lowest value.
To get the series markets risk, as measured by conditional volatility, we model the logarithmic returns of indexes using GARCH models (Bollerslev) 1986). We use the latter econometric technique to measure the systemic risk of equity markets while capturing stylized facts of equity returns, such as volatility clustering (i.c., large (small) swings in stock prices tend to group), hetcroscodasstiitty (i.e., time-dependent variance) and fat-taileduess (i.e., kurlosis greater than 3). With regards to GARCH parameters, select the best combination according to the lowest value of BIC criterion (Schwarz 1978
2hhe website is reachable at httpss//sttoqq.p
TTe risk-adjjsted return is performance metric, defined average compounded return volatility ratio. 4Sortin ratio evaluates risk-adjusted performance financial instrument discounting for its downside standard deviation.
36

---

Published in Transactions on Machine Learning Research (10/2023)
AOR
HSI
NKX
SHC
ST1
TWSE
DAX
FMLB
11EX
CAC
UKX
RTS
BVP
NDQQ
TSX
Avy, Comp. (%)
2.12
15.79
-9.33
-7.43
2.15
31.41
7.76
Volatility (%)
25.94
22.32
23.65
19.26
21.07
19.74
29.25
31.12
30..6
28.79
26.24
12
39. 69
32.46
29.00
Risk-ad].. Ret. (%))
0.14
0.10
0.67
0.44
-0.06
1.43
0.34
0.04
-0.31
0.07
-0.28
-0. 16
0.05
0.97
0.27
Sortino (%6
0.13
0.80
-0.41
0.35
MDD (%6
37.09
25.33
31.27
14.62
31.93
72
38.78
41.54
30.43
38.56
34.93
49.46
46.82
30.1
37.43
ACR/MDD
0.08
0.57
-0.04
0.26
0.03
-0.24
-0.11
1.04
0:21
Skew
-1.10
-0.37
0.27
-0.76
-0.41
-0.54
-0.63
-2.26
-1.05
-0.96
-0.80
-1.02
-1.01
-0.69
-1.01
Kurtosls
7.37
1.87
5.01
6.49
7.22
5.40
10.27
20.19
11.14
9.49
8.94
7.01
11.32
7.35
18:77
%-lle
-6.16
-9.37
5th S-dle (%)
-2.33
-2.28
-2 16
-1.86
-1.75
-1.80
+3.38
-2.81
-2.95
-3.00
-2.79
-3.35
-3.21
-3.12
-2.12
-9.52
-5.56
-6.0 06
-7.72
-7.35
-5.83
112.2
-16.0 92
-14.00
12.28
10.88
-13.02
-14.75
-12.32
-12.34
Max
6.56
5.05
8.01
5.71
6.07
6.37
10.98
8.93
8.57
8.39
9.05
9.23
13.91
9.35
11.96
Table 3: Summary statistics of equity markets at daily frequency, Average compounded return, volatility, risk-adjusted return, and Sortino Ratio are annualised.
37

---

Published in Transactions on Machine Learning Research (10/2023)
Appendix
Additional Results Concerning the Causal Analysis of the Risk of Global Equity Markets
Single scale
Single scale
Single scale
AOR...
AGR NIXK 5TL
AOR... NEX
NAX
STI.
5Th
N15
CAC
FMR
FMIB.
FMIB
UKX
UIC
UKK
NOO
NDO
NDO
AOR
ADR
MKX
NIX
SHC
SHC
SHC.
TWSE.
TWSE
TWSE.
CAG
CAC
CAC
FMIL
FMILL
UIO
UKK
UKK
EVP.
BVP
8VP
TSK
TSX.
T5X
(a) c = 0.01
(b) C = 0.05
(c) è = 0.1
Single scale
Single scale
Single scale
TSX
T5X
TSX
NDQ
ND0
NDQ
8VP
BVP
BVP
UKX
UKX
UKX
IREX
IREX
IBEX
FMIB
FMIB
FMBB
DAX
BAX
DAX
CAC
CAC
CAC
RTS
RTS
RTS
TWSE
TWSE
TWSE
STI
SHC
SHC
SHC
HSI
(d) ē == 0.01
(e) ē = 0.05
(f) 2= 0.1
Figure 17: Highly persistent causal matrices (top) and corresponding SSCGs (bottom) for three different
values of č. Red entries in the causal matrices represent positive coefficients, whereas blue ones are negative. With regards to the SSGS, green nodes are American slock indexes, pink the European ones, and purple the Asian.
38

---

Published in Transactions on Machine Learning Research (10/2023)
(a) ē 0.01
(b) č - 0.05
Scale
5
2
(c) d= 0.1
Figure 18: Highly persistent. multiscale causal matrix for three different values of č. Red entries in the causal matrices represent positive coefficients, whereas blue ones negative.
(a) ē 0.01
(b) 7= 0.05
Scale
Stak
Scat
State
(c) c= 0.1
Figure 19: Highly persistent MSCG for three different values of ē. With regards to the SSGS, green nodes are American stock indexes, pink the European ones, and purple the Asian.
39


Reviews:
[{"review_id": "zjXbuXghcg", "reviewer": "Reviewer_wVxo", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths: \n - Very well-written paper that both identifies a useful generalization of causal discovery for times series (multi-scale data), and gives careful treatment to optimizing the required objective. This leads to an algorithm that is both efficient and effective at recovering causal graphs on simulated data.\n\nWeaknesses:\n - Given its centrality to the story, I would have liked to see a little more background on wavelets and their interpretation in the main text to make the paper self-contained. In particular, **[requested change 1]** I think you should use the finance example to explain why we might expect multi-scale causes (i.e. give some simple examples of causes we might expect for each scale), and then what we'd need to assume for the wavelet transform to be the right tool to detect these potential causes.\n - The experimental section is generally nice, but misses the opportunity to use synthetic data to validate key practical questions. In particular:\n1. **[Requested change 2]** What is the empirical effect of the choice of $\\lambda$ on the performance of MS-CASTLE? You could illustrate this with a plot similar to figure 4.\n2. **[Requested change 3]** How does the performance of the $\\lambda$ selected according to the procedure outlined in section 5.1 compare to the \"optimal\" $\\lambda$ that was selected via hyperparameter tuning (given that in practice you don't have any way of evaluating SHD on real data, this is obviously very important).\n3. **[Requested change 4]** In the real data section, I would have liked to see some discussion on how you selected the wavelet & filter length. Given the sensitivity of the algorithm to these choices illustrated in figure 4, it would be helpful for you to articulate explicitly what you were considering when choosing these parameters. Tying these considerations back to the example used in **[requested change 1]** would practitioners when considering whether to use your method. Given that you don't have an explicit identifiability result, it is worth making these considerations explicit.\n - I would have liked to see a brief discussion on why the discovered relationships are plausible. In particular, it's not obvious to me why volatility in the Canadian market should cause increased volatility elsewhere? How do we know that these markets aren't just experiencing the effect of some common cause? E.g. Brazil & Canada are on roughly the same time zone with roughly the same sized economies, so it may be that effect of global uncertainty are reflected on these exchanges at a different time to the rest of the world. Obviously a detailed discussion of this is beyond the scope of a machine learning journal, but a short paragraph on why it's plausible (or even why it's not!), would make it clearer how these tools should be used.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper introduces an algrorithm that they title, Multiscale-Causal Structure Learning (MS-CASTLE), for causal structure learning in time series systems with multi-scale causal dependencies. The method leverages wavelet transforms to extract causal relationships at multiple time scales, enabling the exploration of causal dependencies across various time lags and resolutions. They derive an efficient ADMM-based optimization routine to solve the problem, and evaluate their method on simulated data. The also show an application of MS-CASTLE to financial market data during the COVID-19 pandemic. The paper also presents a single-scale version of the algorithm, SS-CASTLE, which outperforms other single-scale baselines while reducing computational costs on their simulated benchmarks. The method does not come with identifiability guarantees and is restricted to linear SCMs, but it offers strong empirical performance and computational benefits relative to the baselines.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "See the requested changes above.", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "MzHXjM2it8", "reviewer": "Reviewer_Ze7N", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\nThe biggest strength of this approach is that it is a fairly simple technique that allows for a more robust look at causal inference.  The model is easy to understand, and the optimization approach works for reasonable sized graphs.  Performance is comparable to more complex models despite the limitation of a linear model; however, the linear model is acting on a decomposition, allowing for some additional properties compared to VAR based approaches.  There are many comparisons in the results, showing reasonable performance.\n\nWeaknesses:\n\nI disagree with some of the claims in the introduction.  For example, the authors note \"Consider for example fMRI data regarding different brain regions of interest (ROIs) collected at timestamps t · ∆t, with ∆t being the time interval between two consecutive samples. The aforementioned limit translates into the inability of studying causal interactions between ROIs over time resolutions coarser than ∆t.\"  I disagree with that assertion.  Many of the commonly used techniques in neural data analysis will analyze causal relationships at different scales.  It is well-known that the VAR-based granger causality has an equivalent spectral representation (e.g., spectral granger causality), and that it will look for relationships at different frequencies, and hence different scales.  Likewise, things like phase lag index and directed transfer function often look at spectral decompositions.  I'd recommend looking at [1] for some of the techniques that are used to look across scales in practice.   In contrast, the contribution here is that this decomposition approach allows one to look at connections between the scales, which are not possible in the aforementioned methods.  Hence, the authors need to be more precise in the intro about the nature of their contributions.\n\nThe details of the wavelet transform are lacking, and I would encourage the authors to discuss what wavelets are appropriate.  How are the wavelets being applied to transform the signal in practice?  Does the transformation maintain the causal (e.g., time $t$ only uses data up to time $t$) structure?  Are you limited to wavelets that have orthogonal structure?  Many researchers looking at brain signals, one of the motivating cases, use Morlet wavelets.  Would those be amenable to your framework?  I'd suggest clarifying what wavelets are appropriate.   On a similar vein, how representative are the results given in Figure 4?  While you are exploring different wavelets, they are all in the same family and highly related.  Does this give an optimistic view of the robustness of the representation?\n\nFor the metrics of F1 and SHD, please provide the mathematical definitions and how they are being compared across models since the size of your matrices do not match the VAR-based methods.  There are multiple ways you could be calculating this, and it would be better to have a clearly defined number.\n\nThe procedure for determining the tuning parameters and sparsity strengths is problematic.  Specifically, the authors note “we fine-tune the sparsity strength λ onto separate data sets generated according to the procedure explained above.”  However, in practice we do not have such datasets with known ground truth.  The authors should note this limitation in practice, and highly an approach that could be used whenever labels are not available.  Secondly, the authors note that they choose the value “according to F1-score and SHD.”  However, these don't necessarily agree.  Please clarify.\n\nOne of the things that was not performed was an analysis of how well this method can forecast.  While not directly related to learning the causal structure, the ability to forecast highlights that the method is properly capturing a realistic model.  Here, this could be estimated by inverting the predictions of the wavelet transform.  This is not required for this manuscript, but is a suggestion to help compare the models.\n\n[1] Bastos, André M., and Jan-Mathijs Schoffelen. \"A tutorial review of functional connectivity analysis methods and their interpretational pitfalls.\" Frontiers in systems neuroscience 9 (2016): 175.\n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This manuscript introduces a wavelet-based approach for multiscale causal inference, and then proposes an optimization algorithm to solve for the parameters in the model under constraints and penalties.  The wavelet deconstruction allows the analysis to look at different scales in a natural framework, and those scales can be used simultaneously to model the graph.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Please rewrite the introduction to clarify existing work and practical approaches, and add more specificity to the contributions.\n\nPlease address the wavelet questions mentioned above.\n\nPlease give how the metrics were calculated.\n\nNote the practical issue with the tuning approach, and give an approach that could be used in practice.", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "Ukw2H9c6ZS", "reviewer": "Reviewer_GD69", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths\n1. The writing is well-organized and clear.\n2. The proposed formulation of multiscale causal model and the corresponding estimation algorithm are interesting and well-motivated.\n3. The experiments are extensive, which consist of both synthetic and real financial data. I applaud the authors' efforts in carrying out extensive analysis of their method on financial data, which seems a natural application of the proposed multiscale causal model.\n\nWeaknesses\n1. Identification guarantee were not established for the proposed algorithm, as discussed in Sec. 6.\n2. The method is limited to linear relationship, as also discussed in Sec. 6.\n3. There is limited novelty of the algorithm itself--the key formulation (6) seems like an extension of DYNOTEARS by stacking multi-scale matrices. The primary novelty seems to be an improved optimization procedure with linearized ADMM.\n- While I understand that correctness may be emphasized over novelty for the current venue, the paper should consider explaining in details the similarities/differences of its formulation (6) over DYNOTEARS.\n- When reducing (6) to single-scale version in SS-CASTLE, is it identical to DYNOTEARS (except for the linearized ADMM part)? If that is the case, the authors should consider discussing such connection.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This manuscript considers the problem of learning multiscale causal structure from time series data. A multiscale linear causal model is first formulated based on stationary wavelet transform. Building upon recent works in differentiable causal structure learning, the authors develop an estimation algorithm based on sparse structure, dagness constraints, and linearized ADMM. Although the proposed algorithm does not come with guarantees, extensive experiments on synthetic and real financial data illustrate the effectiveness of the algorithm.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1. Since identification guarantee was not established, the authors could consider at least mentioning some related assumptions that may be needed. E.g., based on my understanding, since least square function is used in Eq. (6), some restrictions on the variances of additive noise may be needed, similar to DYNOTEARS. See https://jmlr.org/papers/v15/loh14a.html.\n2. Related to the first point, recent studies have raised some concerns about the type of optimization formulation used in the proposed algorithm. For instance, will the method be prone to data scaling/standardization (https://arxiv.org/abs/2102.13647, https://link.springer.com/article/10.1007/s11063-021-10694-5) and non-convexity (https://arxiv.org/abs/2304.02146)? As far as I know, the former can be especially important because data standardization appears to be a common pre-processing step in financial data considered in the paper.\n3. Related to the similarity between SS-CASTLE and DYNOTEARS discussed in \"Weaknesses\", does the improvement of SS-CASTLE come from the linearized ADMM, since DYNOTEARS simply uses a vanilla augmented Lagrangian?\n4. The authors could consider comparing to  constraint-based method like PCMCI (https://arxiv.org/abs/1702.07007) in the experiments.\n\nMinor clarification: Below Eq. (3), is $\\epsilon$ supposed to be \"endogenous\" or \"exogenous\" noise?", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This manuscript considers the problem of causal structure learning, whereby we seek to uncover underlying causal relationships in data. This is an important and well-studied problem. The authors however introduce a novel twist on the problem: learning causal relationships that can happen at different time scales. The authors provide a solid mathematical formulation of this general problem and introduce an algorithm to solve it: MS-CASTLE. This algorithm performs admirably in a series of well-designed experiments.

The reviewers universally agree that the work presented in this manuscript is high quality and of interest to the TMLR audience. They are universal in their recommendation of acceptance.

Throughout the review and discussion period, the reviewers provided feedback and suggestions that the authors have faithfully incorporated into the manuscript. I believe these changes have considerably strengthened the work. I commend both authors and reviewers for their engagement throughout this process.

In preparing the camera ready version of their manuscript, I would like to kindly request that the authors take a careful pass over their references, which contain some errors (e.g., "Journal of computational neuroscience," "The Journal of finance," "Cambridge university press," the reference on line 800 is repeated (the first time with "Schlkopf"), etc.).
Decision: Accept (poster)