Learning Multiscale Non-stationary Causal Structures
Authors: Gabriele D'Acunto
Gianmarco De Francisci Morales
Paolo Bajardi
Francesco Bonchi
Venue: Accepted by TMLR ()

Abstract: This paper addresses a gap in the current state of the art by providing a solution for modeling causal relationships that evolve over time and occur at different time scales. Specifically, we introduce the multiscale non-stationary directed acyclic graph (MN-DAG), a framework for modeling multivariate time series data. Our contribution is twofold. Firstly, we expose a probabilistic generative model by leveraging results from spectral and causality theories. Our model allows sampling an MN-DAG according to user-specified priors on the time-dependence and multiscale properties of the causal graph. Secondly, we devise a Bayesian method named Multiscale Non-stationary Causal Structure Learner (MN-CASTLE) that uses stochastic variational inference to estimate MN-DAGs. The method also exploits information from the local partial correlation between time series over different time resolutions. The data generated from an MN-DAG reproduces well-known features of time series in different domains, such as volatility clustering and serial correlation. Additionally, we show the superior performance of MN-CASTLE on synthetic data with different multiscale and non-stationary properties compared to baseline models. Finally, we apply MN-CASTLE to identify the drivers of the natural gas prices in the US market. Causal relationships have strengthened during the COVID-19 outbreak and the Russian invasion of Ukraine, a fact that baseline methods fail to capture. MN-CASTLE identifies the causal impact of critical economic drivers on natural gas prices, such as seasonal factors, economic uncertainty, oil prices, and gas storage deviations.

Published in Transactions on Machine Learning Research (10/2023)
Learning Multiscale Non-stationary Causal Structures
Gabriele D'Acunto DIAG, Sapienza University of Rome Centai Institute, Turin, Italy
gabrisle.dccamloounnhnnnn..
Gianmarco De Francisci Morales Centai Instilute, Turin, Itly
gdff Seffmeeentaiic.m
Paolo Bajardi
Centai Institute, Turin, Iaaly
Francesco Bonchi
Trancesco.oocchi@ceente
Centai Institute, Turin, Italy
Reviewed OpenReview:
https: https:///penreeviau openreview. net/ forum? id= SQmPE633jtA
Abstract
This paper addresses gap in the current state of the art by providing a solution for model- ing causal relationships that evolve over time and occur at different time scales. Specifically, introduce the multiscale non-stationary directed acyclic graph (MN-DAG), framework for modeling multivariate time series data. Our contribution is twofold. Firstly, expose a probabilistic generative model by leveraging results from spectral and causality theories. Our model allows sampling an MN-DAG according to user-specified priors on the time- dependence and multiscale properties of the causal graph. Secondly, we devise a Bayesian method named Mulliscale Non-slalionary Causal Structure Learner (MN-CASTLE) that uses stochastic voriational inference estimate MN-DAGs. The method also exploits infor- mation from the local partial correlation between time series over different time resolutions. The data generated from an MN-DAG reproduces well-known features of time series in dif- ferent domains, such as volatility clustering and serial correlation. Additionally, we show the superior performance of MN-CASTLE on synthetic data with different multiscale and non-stationary properties compared to baseline models. Finally, we apply MN-CASTLE to identify the drivers of the natural gas prices in the US market. Causal relationships have strengthened during the COVID-19 outbreak and the Russian invasion of Ukraine, a fact that bascline methods fail to capture. MN-CASTLE identifies the causal impact of critical economic drivers natural gas prices, such seasonal factors, economic uncertainty, oil prices, and gas storage deviations.
1 Introduction
A causal graph describes causal relationships among the constituents of a given system, and represents powerful tool to analyze such system under interventions and distribution changes. In general, causal graphs are unknown. Fortunalely, it is possible to leverage causal structure learning approaches to unveil and quantify the causal relationships among variables. While randomized experiments are the gold standard for testing causal hypotheses (especially medicine and the social sciencos), in many cases such interventional approaches are unfeasible or unethical. Hence, great effort has been devoted to the development of methods able to retrieve causal structures from observational data [Glymour al. 2019 Schölkopf al. 2021

---

Published in Transactions on Machine Learning Research (10/2023)
Regardless of the different causal structure learning methods, the most informative causal graph is directed acyclic graph (DAG), where the nodes in V are the variables of the system, all possible edges eij E & C VVV are directed and represent direct causal effects, and feedback loops among nodes are forbidden (acyclicity roquirement). A DAG can be associated with its functional representation, also known as structural equation model (SEM, Pearl 2009). Here each node of the causal graph is written as a function of the values of a set of parents nodes and of an endogenous latent noise (see Appendix Ap. In this work, we focus on the case in which such functions are linear and the latent noise is additive.
Even though widely studied and applied, linear SEM is not adequate to cope with causal relations that evolve over time and occur at different time scales, which are both common when dealing with time scries. Indeed, a SEM assumes that (i) causal edges and their weights are stationary and (ii) there is only one time scale at which causal relations occur, ie., the one associated with the frequency of observed data. However, in practice, causal structures might be non-stationary (Zhang et al. 2021b Raggad] 2021 D' AcAuuto et al. 2021 and often there no prior knowlodge about the temporal resolutions at which causal relations occur (Besserve
et al. 2010 Gong et al., 2015 Runge el. al. 2019 D' Acunto el al. 2022).
To overcome these limits, we introduce multiscale non-stationary causal structures, namely MN-DAGs, that generalize linear DAGs to the time-frequency domain. In our work, the term multiscale means that we consider multiple time resolutions, i.e., frequency bands. Hence, we look for causal interactions among time series within each of those distinct frequency bands, and we simultaneously inspect the behaviour of these causal relationships along time. Throughout the paper, we use 2 to represent a certain temporal resolution, where j {L.....J })} indicates the associated scale level and J E N is the maximum level considered. To clarify the meaning of time scale, let consider a data set X E R made by N time series of length T = 2/. The column X[[]= = [X. [f]] represents a sample collected at frequency At. For example, let us say that the time series are observed at daily frequency, hence At is equal to one day. The scale level j = 1 refers to the variations of the time series associated with the time scale of 2 At, i.e., two consecutive days. Analogously, the scale lovel j = 2 refers to the variations of the time scrics associated with the temporal resolution of 2 A1, i.e., four consecutive days. And so on and so forth, until we reach the maximum level J. Additionally, the j-th time scale corresponds to the frequency band [1/2i+1, 11/2/+111/2/).
In MN-DAGs each time scale is represented by different graph page (akin to multi-layer networks). Then, the vertices within a certain page are associated with the multiscale representation of the N time series a the frequency band corresponding to that page. There exists a unique global causal ordering '<' shared by all graph pages, such that the possible parent set P.,- for the i-th node X, can include only those nodes X; that precede it in the causal ordering (X, < X.). Causal relationships among nodes, represented as directed edges, can vary smoothly over time and constitute acyelic structures within each time scale. So, throughout the paper, the term non-stationarity associated with causal structures refers to a smooth dependence on time, similarly to how it is defined by Huang et al. 2020)..
To achieve our goal, the main technical challenges we face are: (i) the definition of a probabilistic generative model that allows to sample an MN-DAG; (ti the development of a learning method to estimate MN-DAGs from real-world data. Regarding the first point, we propose a probabilistic generative model over MN-DAGs having the causal ordering and the causal relationships as latent variables. The observables are N zero- mcan time scries of length T. Our model leverages lincar SEM and multivariate locally stationary processes (MLSW, Park et al. 2014). In particular, MLSW is mathematical framework to represent time series as sum of contributions coming from different time scales (see Appendix C). Concerning the second point, we propose a Bayesian method, called MN-CASTLE, that uses stochastic variational inference (SVI, see Appendix D for learning causal structures. MN-CASTLE is able to cope with multiscale data that features time-dependent variance. Our method relies upon observational data and the estimate of the inverse of the power spectrum at different temporal resolutions.
Overall, our contributions can be summarized as follows:
We define new type of causal structure for modeling causal relationships that evolve over time and occur at different time scales (MN-DAG).
2

---

Published in Transactions on Machine Learning Research (10/2023)
We devise a probabilistic generative model that allows sampling an MN-DAG according to user-specified priors on the time-dependence and multiscale properties of the domain. Our model can be used to generate synthetic time series with real-world characteristics.
We design a Bayesian inference method, MN-CASTLE, for estimating MN-DAGs from real-world data. Our empirical assessment on synthetic datasets demonstrates that MN-CASTLE outperforms baseline methods in various experimental settings and is robust to model misspecification.
When applied to study what drives natural gas prices in the US market, MN-CASTLE succeeds where baselines fail. In fact, it is the only method able to capture the dynamic nature of the market and the impact of exogenous events such as COVID-19 and the Russian invasion of Ukraine.
At a high level, this paper bridges the gap between multiscale modeling and machine learning- based causal structure learning methods.
Roadmap. This article is organized follows. [eection 2 relates our proposals to existing methods, highlight- ing differences and similarities. Then, the subsequent three sections deal with our first technical challenge. Specifically, Secciio 3 presents our probabilistic generative model, [Section 4 details how an MN-DAG is sampled in our model, and [ecctio 5) how to generate data from the sampled MN-DAG. At this point, we address our second challenge in [Section 6 where we pose a Bayesian learning method developed according to the proposed probabilistic generative model. Next, [Section 7 presents the empirical assessment of our model. In detail, [Section 7.1 statistically describes data generated by the probabílistic generative model. Section 7.2 regards lests on synthetic datasets, by providing details concerning the experimental settings and introducing the considered baseline models. Subsequently, [Section & analyses a real-world use case natural gas prices in the US market. Finally, [Section 9 concludes with an additional discussion concerning our findings, and outlines open questions and future research directions.
2 Related Work
Causal structure learning methods can be mainly classified into three categories, according to the approach used to infer the causal graph: (i) constraint- approaches, which make use of conditional independence
tests to establish the presence of a link between two variables (Spirtes et al. 2000 Huang et al. 2020 (ii)
score-based methods, which use search procedures in order to optimize a certain score function (Heckerman)
et al. 1995 Chickering 2002 Huang et al., 2018|; (ifi) functional causal models, which express variable
certain node as a function of its parents (Shimizu et al.) 2006 Hoyer et al., 2008 Zhang & Hyvärinen 2010;
Shimizu et al., 2011: Peters et al., [2014] [Bhhlmann et al. 2014). Our approach fits into the latter calegory and aims to handle the presence of non-stationarity and different temporal resolutions in the underlying
causal structure.
Unlike the multiscale causal structure learning method proposed by D' D'Acunto el al. (2022). which estimates multiscale stationary causal relationships hinging on stationary wavelet transform (Nason & Silverman 1995) and non-convex optimization, our method applics a different learning scheme and is able to handle non-stationary relationships
well. Furthermore, the method we propose exploits the estimate of the decomposition of the inverse power spectrum at different time scales, whereas the algorithm proposed in the previous paper operates on the estimated wavelet detail vectors.
In the past, several approaches have been developed that can infer causal structures in the presence of non- stationarity under certain assumptions (Song et al. 2009 Ghassami et al. 2018 Strobl 2019 Perry et al. 2022). The main (implicit) assumption common to these approaches, concerns the time scale at which causal interactions occur, that is, it is assumed that this scale coincides with the frequency of observation of the data. The model we propose, relaxes this assumption, and allows time-dependent causal relationships to be investigated at different temporal resolutions. Another difference concerns the assumption regarding the existence of nultiple domains, where causal dependencies between variables may vary but are assumed to be stationary within each domain, to exploit zon-stationarity and distributional shifts to recover the underlying causal structure. Although in the context of time scries, the datasct can be segmented into different domains through a sliding window approach, this procedure introduces discretionary choices such as (i) the choice of the splitting points and (ii) the size of the time window in which causal relationships should be stationary.
3

---

Published in Transactions on Machine Learning Research (10/2023)
However, in general, for real data there is prior knowledge regarding the above issues: the causal structure might vary a lot even when windows are overlapping D'Acunto et al. 2021). In contrast, our method aims to learn the causal structure and describe its temporal evolution, assuming that it is linear in the frequency domain and that the causal ordering is shared between the temporal resolutions considered.
Our probabilistic generative model extends the works of Cundy et al. (2021;: :: Charpentier et al. (2022) , since it is suitable for time-series data and provides a causal structure that lives in the time-frequency domain. Even though our approach leverage Gumbel distributed variables for sampling the causal ordering in the previous two works, the procedure we apply is different and requires lower computational cost Gadetsky et al. 2020). In addition, our inference model uses gradient estimator with a data-dependent control variate strategy for learning the parameters of the causal ordering distribution, whereas existing models exploit differentiable relaxations of such a distribution. Our procedure uses the masking of distributions as well to optimize at each step only the causal relations compliant to a certain causal ordering. A similar approach is also employed in Ke et al. (2019); Ng et al. (2022). However the masking used in those works aims at excluding all non-causal relations, not just. those that do not conform to the causal ordering.
Finally, exploit recent developments in variational inference in order to approximate the posterior dis- tribution over MN-DAG parameters given data, in accordance with the MN-CASTLE probabilistic model. This general learning scheme is also exploited in other recent works (Cundy et al., 2021 Charpentier et al. 2022 Aunadani et al. 2021 Lorch et al., 2022) to model the posterior distribution over the parameters of a DAG, as defined in the correspouding proposed probabilistic models.
3 Probabilistic Generative Model
The probabilistic generative model over MN-DAGs we put forth incorporates both the causal ordering and the causal relationships as latent variables. The observables consists of N zero-mean time series, each of length T. In our model, the underlying causal structure determines the transfer function matris. Specifically, this matrix is defined as a time-dependent mixing matrix that depends on the causal structure and the strength of causal relations a measure of the local variance and cross-covariance between the time series, which is equivalent to the power spectral matrix (see [Section 5). From a lincar conditional dependence perspective, the inverse of the power spectral density is a concept extensively studied in spectral theory (Dahlhaus 2000 and references therein). Indeed, the inverse of the power spectral density is generalization of the precision matrix to the time-frequency domain. It provides information about the local linear dependence between two time series after removing the linear effects of the rest of the time series. Therefore, according to our model both the power spectral matrix and its inverse are driven by the multiscale time-dependent causal structure. This dependence on the causal structure implies that the power spectral matrix and its inverse are also time-dependent.
each time step and scale level. Furthermore, the transfer function matrix provides
The proposed probabilistic model takes as input (i) the number of nodes N E N; (ii) the number of samples T E N; (iii) a parameter f E [0, associated with the multiscale feature; (iv) r E [0, I] which describes the time dependence of causal relationships; (v) 8 E [0, I] that manages the density of the MN-DAG.
Figure shows a (T, p)-quadrantt along with examples of latent causal structures which determine the sampled data, according to the specified values of μ and 5. When μ = 0, we obtain the single-scale case depicted in Figure 1[a) and [bb. Here, the MN-DAG has only one page. We assume that the power spectrum that describes the system is concentrated in the finest scale level j = 1. In addition, if also T 0, the causal links are stationary Figure a()).. Starting from the origin, as we move to the right (T + 1) the temporal dependence of causal connections increases. As we move upwards (u + 1), the likelihood that the causal graph contains more pages increases FFiure T(c) and dd)) Then, the overall power spectrum is spread over more temporal resolutions.
The following Sections and 5 delve into the sampling of the MN-DAG and the generation of data, respec- tively.
4

---

Published in Transactions on Machine Learning Research (10/2023)
time-scales
[1,0]
(1.1)
(multiscalo)
(0,0)
[0.1)
(nonstationarity)
2
b)
a)
Figure 1: Examples of causal structures that can sampled from the proposed probabilistic model according the specified values for not-stationarity and multiscale features. In the depicted quadrant, we have the non-stationarity feature (associated with the parameter feature (associated with the parameter μ [0, 1)) the y-axis in yellow, respectively. Colors, edges shape and number of graph layers highlight differences from the single-scale stationary DAG corresponding to the origin of the quadrant (a). When we move horizontally, the temporal dependence of the causal coefficients (odges in the cansal graph) changes (b). Similarly, vertical shifts in the quadrant are associated to the change in the number of time scales (pages the causal graph) contributing to the sampled data (c). Finally, when both 7 and μ are different from zero, we sample data concerning a system driven by an underlying multiscale non-stationary causal structure (MN-DAG).
[0, 1]) the x-axis in red, and the multiscale
4 Sampling an MN-DAG
Figure 2 shows the three steps needed to sample an MN-DAG.
Sample the time scales. The number of pages (time scales) of the MN-DAG is 1 J', where J is sampled from binomial distribution J' ~ ((log2(T) 1). Here, the first parameter of the binomial distribution is the number of trials and represents the probability of success. Without loss of generality, we assume that temporal resolutions are consecutive, i..,, given the value of J, all the time scales 2/, JI, are associated with a page in the causal graph. This assumption does not imply that causal relations occur within all the considered pages. Since the model probabilistic and the user specifies a value for the density 8, we also might end up with a causal graph without edges.
Sample the causal ordering. Within our probabilistic model, we assume that the causal ordering is shared by all time scales. This property implies that, given x, the possible parent sets at each temporal resolutions P.x for the i-th variable X; (P: P. < X.). The causal ordering < be thought
5

---

Published in Transactions on Machine Learning Research (10/2023)
a) Sample time-scales (pages graph)
b) Sample the causal ordering
MN-DAG
Causal ordering < (2((8)
= 1 J J' ~ B (T)
j
1, 1b)
13,
j=1
j+3
3
j=1
Sample the causal tensor
j=2
c.1) Sample weight tensor
c.2) Sample logical
c.3 Causal
c.4) Impase the causal ardering
tensor II
coofficients tensor C
Suue
[(6)
C P'PP
3
j+2
j=3
J-3
3
Figure The figure shows the steps necessary sample a MN-DAG. For the sake of readability, let us consider the case where N-3 and T-4. Here the yellow color refers time scales; red for non-stationarity. (a) First, sample the number of pages (time scales) of the MN-DAG. Given μ4, the latter is given by J= J', where J' ~ B(Tog2(T) - 1,)... In the example, we instantiate three pages (J 3). (b) Second, we sample the causal ordering <~ PL(@) that is shared by all time scales and entails the permutation tensor P (0, (JJJTXXNN Here, PL(0) indicates the Plackett-Luce distribution, defined by a score vector E RN,, where 0, ~ N), The indexes are j for time scales, t for time steps, n for the considered nodes and n' for the positions within x. the considered example, we have <= [3, 1]'. Therefore, for each 3 slice of the tensor corresponding to certain time scale j and time t, we have Pn'n = 1 (blue square) the node appears at index n' within <. (c) Third, build the tensor of causal coefficients follows. With regards indexes, j, & and are the same as above, whereas indicates the parents dimension. Here, we first sample a full tensor of weights W E J made by three components: (i) a constant term WD; (ii) W that makes the magnitude of causal relationships different across time scales; (iii) W that allows causal coefficients to vary over time according to batched GP(O,K), Therefore, within each scale j, Wam manage the density of the entailed MN-DAG, we multiply
are smooth functions varying over index element-wise W by a logical mask II {0, 1 The entries of the II are distributed according Bernoulli distribution, 7mm B()). Finally, obtain the causal tensor C that entails the MN-DAG the right by imposing the causal ordering sampled at step (b).
permutation of a vector of integers <'= [I,...,N],
N). thus we use the Plackett- Luce distribution (PI,, Plackett 1975 Luce 1959 to sample it. PL represents a distribution over permutations, defined by a vector scores 0 R which allows sampling permutations Sx is the support of permutations of N elements. Thus, given 0, the probability of permutation b is
Sy in O log N), where vector of N integers and
p(b 0)
II
A sample from PL distribution can be thought as sequence of samples from categorical distributions: first by comes from the categorical distribution with logits 0; by from the categorical with logits 0 - {0b, and so on. The mode of the PL is the descending order permutation of scores bo = 0gg 2 Agg 2 2 0gg The sampling procedure from PL relies upon the fact that order of vector R Gum((.1) is

---

Published in Transactions on Machine Learning Research (10/2023)
distributed as PL(O), where Gum 1) is Gumbel distribution with location parameter @ and scale equal to one (Gadetsky et al. 2020). Therefore, we can sample b as follows:
2, 0, log( log(e), (z)) = argsort(z).
2 7~~ (0,1)
We sample the causal ordering Y~ PL(O) by using the procedure above, where choose a uniform prior for the PI score vector, i.e., O₂ N). The causal ordering < entails permutation matrix P (0, NNNNN such that Pa'n iff the variable X,, occurs at position n' within <; 0 otherwisc. Finally, derive permutation tensor P E R X by simply tiling P along both multiscale and lime dimensions.
Sample the causal tensor. Given J and T, we build a lensor of weights W E RJXXXXXN made by three building blocks. First, we sample W E RNxN whose entries are normally distributed, wh ~ (N0 1). Starting from W we derive the first component w E / XNNN by simply expanding W along both multiscale and time dimensions. Then, this component can be thought as a constant term shared by all temporal resolutions and time steps.
Second, we sample W E IxIxNNN whose entries are distributed according to a Gaussian (0,, By ex- panding T along the time dimension, we obtain the second component W E IR /XTXNXN This component makes the magnitude of causal relationships different across scales and is stationary along L.
Third, we sample W E JxTNNXN where each tube along the time dimension follows a multivariate Gaussian distribution MN(O,K) K). Here, the covariance matrix valid kernel(s) for Gaussian processes (GP, Bishop & Nasrabadi 2006), where the lengthscale is À 1/T. This component imposes the causal coefficient to evolve smoothly over time, according to Indeed, T - 0, the lengthscale of the kernel increases and consoquently W varies more slowly in time. Finally, the tensor of weights is
K(t, represents (combination of)
W = w + W W
(1)
Now, to manage sparsity and ensure the acyclicity of cansal connections, we generate a suitable logical mask EE (0, (0,11)a 1JJXIXNxN, Indeed, we use this mask to obtain a tensor of causal relations from W, made of strictly lower triangular slices. The slices Nam are strictly lower triangular and the entries are distributed according to a Bernoulli distribution, Tum ~ B6)). Then, we obtain the tensor of causal relations as Č II o W, whose slices Crmm are nilpotent Here II {0. 1X.XXXXXN is obtained by expanding II over time and o represents the Hadamard product.
At this point, given P and Č, we compute the causal tensor that entails the latent MN-DAG by means of the product C P'CP, where P' is obtained by transposing the two rightmost dimensions of P.
5 Generate Data from the MN-DAG
Having sampled an MN-DAG, we wish to use it to generate N zero-mean processes of length T, whose behaviour is determined by the evolution over time of a latent MN-DAG. Here, we build upon the SEM App- pendix A) and the MLSW (Appendix C theoretical frameworks. Mathematieally, we model the multivariate time series as
+oc
x
Σ
Σ M,lu8,,,,,3 k].
(2)
In Equation (2) (i) (1 - k] is a set of non-decimated wavelets (see Appendix Bb: (ii) {2,,0} a set C, of random vectors Z,,k INx.N: (iii) M, (I-C, (pir) is time-dependent mixing matrix that represents the transfer function matrix, where = k/T is the rescaled time (Dahlhaus 1997) and
A matrix A is said nilpotent if is square and A integers N, where N known the index of A.

---

Published in Transactions on Machine Learning Research (10/2023)
Cy[u] E RR-NN is the matrix of causal coefficients at time and scale described in Section 4 model, the local variance and cross-covariance between the processes at certain time y and scale J, i.e., the local wavelet spectrol matriz (LWSM) S,[²] is determined by the MN-DAG structure. The same holds for the inverse of the LWSM, i.., Oylv] = S, that provides information concerning the local linear dependence between two time series after having removed the linear effects of the rest of the time series. Thus, we can think of it as the equivalent of the precision matrix the time-froquency domain. Then, the matrix Oz[r] relates to partial correlations between time-series, concept that has been extensively studied in spectral theory (Dahlhaus 2000 and references therein). Indeed, rescaling leads to the partial coherence
In our
R = D)I/oo,\viDDbb) with Dy(v] being a diagonal matrix with entries [Oy[v]] Yn E [M] Given delinition, the partial coherence between two time series provides measure of local linear dependence as well and is bounded in [-1,1]. The results below link the spectral properties of XT to the underlying multiscale causal structure.
Lemma 5.1. The transfer function matrix permuted lower triangular matria, M, P'(C-, where P E NNN a permutation matriz such that Pw'n = 1 iff the node Xn occurs at position n' within the causal ordering x, and C,(v] strictly lower triangular matrie of causal coefficients.
Proof. See Appendix
Equipped with Lemma 5.1
we obtain the following.
Lemma 5.2. The local wavelet spectral matriz and its inverse are given by S,[M] = P'(L - ((U--,,i)~-(
C/ and 0,lo] = P - CI
Proof. See Appendix
Lemma 5.2 provides with the expressions of Sy [r] and Oy[r] in terms of the features of the MN-DAG, that are the causal ordering eutailing the permutation matrix P and the matrix of causal coefficients. From computational perspective, the expression O,l is more appealing than that of S,(²] since it does not involve any matrix inversion. Therefore in Section 6 we leverage Oj[v], given its convenient form and the information that it provides.
It is interesting to understand how the causal ordering and the order in which we observe the dimensions of XTr impact the spectral properties of the process, i.c., the information contained in S,[²] and O, [v]. Proposition 5.3 shows that any permutation of the causal matrix leaves the spectral properties unchanged. In addition, it proves that any re-ordering of the dimensions of XT results in a representation of the form in Equation (2) with the same spectral properties.
Proposition 5.3. The spectral properties of the process XT are independent of both the causal ordering and the order in which the process dimensions are observed.
Proof. See Appendix
6 Two-Step Inference
We expose a Bayesian method for the estimation of MN-DAGs, termed MN-CASTLE. It is implemented by using Pyro (Bingham et al.) 2019: a probabilistic programming language built. on Python and Py- Torch (asszke et al., 2019). A probabilistic model is a stochastic function that generates data x according to latent randoi variables z and parameters 3*, having as joint density function
P-- (x, z) = P3- (x z)pa- (z).
where pa- (z) and Pa- (x z) are the prior and the likelihood, respectively. The goal is to learn the parameters of the model 3 from data. As detailed in Apppendix DJ SVI offers scheme learn 3 by approximating the usually intractable posterior distribution pa- (z x) by means of tractable family of variational distributions (z),, called guides, parameterized by the variational parameters Φb
8

---

Published in Transactions on Machine Learning Research (10/2023)
Our task is as follows. We are given dataset X' = (Xr(t))<1 Xy(t] [x. [t]]]' and estimate of the inverse LWSM O, different time scales j. As an example, the inverse smoothed bias-corrected raw wavelet periodogram is suitable non-parametric estimator (Park et al. 2014). Then, according to the probabilistic generative model in [Section 5] want to learn the following parameters given previous inputs by means of SVI: (i) the vector of scores 0 of the Plackett-Luce distribution used to model the latent global causal ordering <; (ii) the mcan and kernel parameters of the latent batched GPs used to model the entries of the hidden causal coefficients tensor C, i.., of GP((C ((t,t')), Here, we assume that functional form of the kernel ((t is shared by all causal the coefficients. In addition, by learning the kernel parameters, obtain estimate f of since we assume T 1/A in Section 4
In light of Lemma 5.2 and Proposition 5.3 we set the inference of the causal ordering apart from that of the causal coefficients. Indeed, our inference procedure is as follows:
Step 1: We estimate the parameter vector of the PL distribution which determines the causal ordering by conditioning on the dataset N;:
Step 2: We estimate the parameters of the kernels associated with the causal coefficients by conditioning on the estimate Oj, while imposing the causal ordering from Step 1.
Plackett-Luce
Plackett-Luce
C Normal
C Normal
Xr Multivariate normal
R
Açe ER
Oco E R,,
index
inclex
X
Index
0
index
Index mn.
(a) Probababilistic model
(b) Guide
Figure 3: Graphical models associated with a the probabilistic model and b the parameterized variational distribution for learning the causal ordering, along with variational parameters and their constraints.
Model and guide for causal ordering inference. Figures 3a and 3b provide pictorial representation of the probabilistic model and the guide used in Step 1. In particular, we resort to graphical models to illustrate the corresponding joint distributions. Here, random variables are represented as circular nodes, where a blank node represents latent variable while a grey one is associated with an observed variable. Deterministic variables are represented as rhomboid nodes, while variational parameters are printed outside of the nodes. Edges indicate dependence among variables and rectangles (plates) indicate conditionally independent dimensions, i.e., independent copies. In addition, Figure 3 provides the distributions (along with the constraints of parameters) of random variables and variational parameters.
Figure 3 shows that model and guide share the same latent variables. Indeed, since the guide is used to approximate the true posterior, it needs to provide a valid probability density over all hidden variables. More in detail, we have two latent variables: (i) the causal ordering, which is global since it does not depend on any other variable and is modcled within the guide as a PL(@); (ii) a stationary single-scale causal structure Co E NNN, where each entry Co is independent. of the others and is modeled the guide with Gaussian. Since we assume the causal ordering (i) shared by all time scales and (ii) stationary; we infer it from observed data X7, without any additional information concerning the variance decomposition and its evolution over different temporal resolutions (provided by given estimate of S,). For this reason, to learn 0, we resort to the SEM formulation given in Equation (4) where we set C PPCPP (see Section 4). As consoquence, the causal tensor C in Figure 3 depends on x. According to the probabilistic model in Figure 3a] at each time-step we observe the vector X by using multivariate normal likelihood,
9

---

Published in Transactions on Machine Learning Research (10/2023)
precisely (O MM') Indeed, with constant causal coefficients and normally distributed noises, we have:
E [ = M. B. [[] M. 0 0; (ii) Var(Xx(4]] Var[MZ] MTM' = MM'.
Co
(u)
K
index
En((n
indlex
index
index
index
Multivariate normal
index j
~ Normal
Ç E R
Mq(u) E E
Engtu) E R
u ~ Multivariate normal
C ER
Normal
E IR
Ô Normal
A R
(a) Prohababilistic model
(b) Guide
Figure 4: Graphical models associated with a the probabilistic model and b the parameterized variational
distribution for learning batched GPs, along with variational parameters and their constraints.
Model and guide for batched GPs inference. Figures 4a and 4b depict the probabilistic model and the guide used in Step 2. Here we exploit the result from Step We compute the mode of the PL distribution x0= argsort(-O) 0), and then build the permutation matrix P described in Section Since known, represent 20 with rhomboid node.
To model each latent causal coefficient at certain temporal resolution level j as a smoothly varving function f f~ GP(C) k(t,t')), exploit a variational formulation of Gaussian processes (Hensman et al. 2015
Accordingly, we consider set of inducing points {(+)f optimised over the training set, where T < T, and latent inducing function variables (a subset of f) over these inducing points. Thus, the method relies upon the introduction a joint variational distribution (f,, u), such that it factorises as p(f u))(u) (Titsias 2009). This allows avoid the computation K within the inference procedure. Here, to approximate the true GP prior p(u) over the inducing points, we choose q(u) to be a Cholesky variational distribution, i.., multivariate normal with positive definite covariance matrix M MN(Hggu)) Eg(u)). This variational approach allows to reduce the computational burden of GP estimation while avoiding overfitting the same time (Bauer et al. 2016). In detail, in our work the usage of T inducing points lowers the compulations ional cost of each GP from ( to (TT)) (Hensman et al. 2015). As a consequence, in both Figures 4a and 4b we have two latent variables: associated with inducing functions and f associated with GP prior values. Here, we use batched GP to model causal coefficients, consequently they are independent both within and among time scales (rectangles in Figure 4). Since the joint distributions within the model and the guide q(f, u) factorise as;
p(f, u) p([ u)p(u); (f, = pff u)g(u),
also draw edges from to f. The variational parameters, along with their constraints, are shown in Figure 4b
Now, let us consider a lower triangular matrix of ones L, having size N x N. Furthermore, define R, E NXN such that the entries [Ryinn' ||| <E RRymn /T///// Hence, mask the distribution of the hidden functions with B, (P'LP) R, to take into account and update only the relationships conforming to -0 and associated with nonzero partial correlation, where represents the Hadamard product. In practice, since it is difficult to exactly estimate zero partial correlation, before constructing the mask it is possible to hard-threshold the values of R, at a certain value [0. 1]. In our experiments we set P 0.05. At this point, we observe the estimated Ö, by using Gaussian likelihood N([I CI, where the scale E R. is fixed (here 0.05). In particular, the mean value of the latter Gaussian is set
C,),
10

---

Published in Transactions on Machine Learning Research (10/2023)
in accordance with Lemma 5.2 To implement these probabilistic model and guide, we combine Pyro and GPyTorch (Gardner et al. 2018 an efficient Python library for GP inference built on Py
Inference. We optimize the variational parameters above by using SVI, and adopting as optimizer Adam (Knngma & Ba 2014) along with learning rate decay and gradient clipping (Goodfellow et al.) 2016). These tricks are useful to avoid bouncing around local optima when you are close to them and to prevent the gradient from becoming too large. In Step 1, we optimize w.r.t. 0 to approximate the likelihood of given XT. Unfortunately, the latent. variable < is non-reparameterizabl.. Therefore, we use the REINFORCE estimator (Williams, 1992.. which is suitable for getting Monte-Cari estimates of certain cost function f (z). According to REINFORCE, we have
VelBenne) = (x) [(V (ogg(z))f((l)) + Vofo(z)].
(3)
Although unbiased, this estimator known to have high variance. A way for roducing this variance is by means of control variate stralegies, i.e. by adding function within the expectation operator in Equation (3) that depends on the chosen values for z but is constant w.r.t. Φ. So, the additional term does not affect the mean of the gradient estimator. Here, The rationale behind the usage of basclines, is to reduce the variance by tracking the mean value of fo Thus, we add running average of (z), namely fa(z), for predicting the value of fø(z) at each step. On the contrary, in Step 2 we exploit the reparameterization trick (see Appendix D), to the benefit of learning, Finally, we return the learned variational parameters once the maximum number of iterations is reached.
resort to data dependent baseline (Muih & Gregor 2014
7 Results
In this section we present the empirical assessment of our proposal. We first dive in the statist ical analysis of the time series generated by the proposed model in (ecciion 7.1 Then, Section 7.2 presents results regarding the inference of MN-DAGs from synthotic data.
7.1 Probabilistic Model over MN-DAGs
We start by illustrating the output of the proposed probabilistic generative model by means of an example. We consider N = 3 nodes (time series), T 512 time steps, multiscale level μ 0.5, non-stationarity level 0.5, and density of causal interactions 8 = 0.5.
First, Figure 5) displays the underlying MN-DAG, sampled as detailed in Section along with the evolution over time of causal relationships. We obtain an MN-DAG composed of three pages, corresponding to temporal resolutions 2/, j = {1,2,3) The sampled causal ordering is <= 3, 2). and all causal relations, here locally periodic functions with increasing variations, are compliant with x. Indced, we can only observe directed edges from time series n to m, where n < m.
Now, given the sampled MN-DAG, we generate data according to Haar wavelet (Nason et al. 2000) oscillatory function 1, k). Figure 6 depicts the generated time series along with descriptive statistics.
Equation (2) where we use non-decimated
On the first row, we have the behaviour over time of synthetic data. Here, we resort to the augmented Dickey- Fuller test (ADF, Dickey & Fuller 1979) to assess stationarity. According to the test, the null hypothesis Ho indicates that the process has unit root (i.e., is non-stationary). The resulting p-values prove that the generated processcs are (wcakly) stationary (we reject Ho). Indeed, they have zero mean, while their dispersion looks different. On the one hand, the variance of X, (which occur at first position in x) is stationary, on the other hand those of X₂ and X3 vary over time. Furthermore, X2, that has incoming causal edges at all temporal resolutions, displays the largest swings.
On the second row we provide the histograms of the data, where we employ a Jarque-Bera test (JB, Jarque & Bera 1987) to assess normality. In particular, the mull hypothesis Ho is that the process is normally distributed. The resulting p-values suggest that X, is normally distributed, while we reject H, for both X₂ and X3. Indeed, the associated distributions are leptokurtic, with X3 having a more pronounced negative fat tail.
11

---

Published in Transactions on Machine Learning Research (10/2023)
Scale
Scale
Scale 3
C
C
C13 31
500
0.0
55
500
500
2.5
500
Figure 5: The figure depicts the latent MN-DAG sampled by using the proposed probabilistic generative model, where we set the number of nodes N 3, number of time sleps T 512, multiscale level μ = 0.5, non-stationarity level r = 0.5, and density parameter 5 = 0.5. The resulting MN-DAG has (i) 3 nodes; (ii) J = 3 pages (yellow rectangles), (iii) non-stationary causal interactions (red directed arrows, values shown time series in the insets) that follow Gaussian process with kernel K = KPeriodic + Kineear X K 3atern/aa (iv) global causal ordering x= [1,3, Within each scale, we also plot the evolution of causal relations over time. Kernel variances are Linenr Periodic Matern3/2 = 1; the lengthscales Aperiodic = Matern3/> = 1/T, and the period PPeriodie = 1/7. Given the kernel shape, the causal coefficients are locally periodic functions with increasing variation.
Looking at the autocorrelation (with lag / E [0] plotted the third row, we see that all the generated time series show serial correlation, statistically significant al 95% level (light blue bands). This result is in accordance with the multiscale nature of the time series. In addition, the autocorrelation is driven by the local wavelet spectral matrix Sj (see Appendix C) that in our model is determined by the causal structure.
Finally, the autocorrelation of absolute values the processes prove that large swings in X2 and X3, either negative or positive, tend to be followed by other large swings. This effect is also known as volatility clustering, key-feature of financial time series (Mandelbrot, 1967 Ding & Granger 1996). Here, large movements in the series are driven by the increase of causal coefficients modulus, shown in Figure 5)
7.2 Causal Structure Learning from Multiscale Data with Time-dependent Variance
We next report a comparison between our method, MN-CASTLE, the algorithm for multiscale causal struc- ture learning introduced by D' Acuuto et al. 2002), MSCASTLE, and state-of- -he-t---tt--rttt algorithms for learning Equation (4) For this comparison we use basclines belonging to different familics, and synthetic data gener- ated by the proposed probabilistic generative model. The goal is to assess the gain, in terms of performance, we deviate from the single-scale stationary case, ie.,, / = 0, which is the closest to Equation (4) Additionally, we report results concerning the inferred causal ordering.
Settings. We rum our experiments according to four main different configurations. First, to evaluate the methods as we move within the (T, qu-uudrant, we generate the data by setting N = 5 and T = 100, while the entries of the PL score vector are drawn from a uniform distribution 8 ~ (O, N).. . We test three values each for the multiscale and non-stationarity parameters, thus giving raise to configurations of none, medium, and high values for each parameter. For each possible combination (T, {0.0, 0.5, 0.9] {0.0, 0.5, 0.9), we generate 20 datasets that contain N time series each of length T. With regards the causal structure density, use 8 0.5.
Second, to measure the sensitivity of the performances w.r.t. network density, we set (N,T, equal to (5. 100, 0.5, 0.5) and let. 8 varies in (0.25, 0.5, 0.75). For each possible combination, we generate 20 datasets.
Third, to measure the sensitivity of the performances w.r.t. network size, we set (T, T, 14 equal (100, 0.5. 0.5, 0.25) and let N varies in (5, 10, 15, 20). Thus, in this experimental context we go from a configuration in which the number of observations T is greater than the number of relationships possible
12

---

Published in Transactions on Machine Learning Research (10/2023)
X p-val of ADF test:0.000
X2. p-val of ADF 1est:0.000
p-val of ADF test:0.
Value
Value
Value
0-
0-
-20-
-10
200
400
200
400
200
400
Histogram X₁, p-val of test:0.313
Histogram of X p-val of JB test:0.000
Histogram of Ng. p-val of JB test:0.000
Skew:0.00 080
Skew:00001
Skew:-0 166
Count
Kurt:-0.280
Count
Kurt:2.99 910
Count
50
Kurt:1.000
25
5.0
2.5
0.0
2.5
5.0
20
20
10
3
Value
Value
Value
Autocorrelation X,
Autocorrolation of X2
Autocorrelation of X3
0.25
0.0
0..2
0.25
10
20
30
40
10
20
30
40
10
20
30
40
Lag
Lag
Lag
Autocorrelation of Xx
Autocorrelation of X2
Autocorrelation of Xa
0.1
0.2
0.25
0.0
0.00
0.0
10
Lag
Lag
Lag
Figure 6: The figure shows the generated time series, along with descriptive statisties, where each process is associated with different column. (i) Starting from the top, we have the synthetic data obeying to the underlying MN-DAG, where we provide the p-values of an ADF test. (ii) On the second row, we have the histograms of observed values, along with the p-values of a JB lest, skewness, and kurtosis. (iii) The third row shows time series autocorrelations (with lag / (1....,,00) 40))- The light blue bands show 95% CTs. (iv) The last row shows the autocorrelations of absolute values of the processes.
13

---

Published in Transactions on Machine Learning Research (10/2023)
in complete single-scale DAG, iee, N (N- 1)/2, to one in which it is much less. Also in this case, we generate 20 datasets for each combination.
Fourth, we look at the performances of MN-CASTLE under generative model misspecification, ie,, when the assumptions underlying Equation (2) are violated. Here, we consider two kinds of violations. First, we set (N,T, T.) equal to (5, 100, 0.5, 0.5) and violate the Gaussianity assumption of the latent noise 2.J.- Indeed, we generate 20 datasets for the case in which žj,k ~ L(00 1) and 20 for that where Zj,k ~ U(, 1), being 1) the Laplace distribution with zero mean and unit scale and (O,,1) 1) the uniform distribution over the unit interval. Second, we run MN-CASTLE on downsampled data, meaning that the frequency associated with the observational task is lower than that at which causal interactions occur. In detail, we generate 20 datasets by setting (N, T, H) equal to (5, 256, 0.5, 0.5). Then, we subsample the observations to T = 128 and run MN-CASTLE on these decimated datasets. This way, our model cannot exploit the information related to the scale level j = 1, and consequently cannot infer the causal structure at that time resolution. Hence, we test whether the lack of such information affects the performances of MN-CASTLE in retrieving the causal relationships at coarser time resolutions.
Finally, in case * 0, for the GP we use the radial basis function kernel KRBF with variance OBBF 0.1 and lengthscale ARBF= = 1/T.
Baselines. We test MN-CASTLE against the following four baseline models. First we consider MSCAS- TLE, a multiscale causal structure learning model which exploits multiresolution analysis and non-convex continuous optimization to retrieve stationary causal relationships. Next, we have DirectLiNGAM (Shimizu jt al. )0 & method belonging to the family of non-Gaussian models. Algorithms within this class assume that the noise Z is non-normally distributed. Indeed, in this case the causal structure has shown to be fully identifiable (Shimizu et al. 2006). Then, DirectLNNAAM returns an estimation of both causal ordering and causal coefficients. Second, we have CD-NOD (Huang et al.) 2020), which belongs to the family of constraint- based methods. In particular, it has been developed to deal with heterogeneous (no assumptions on data distributions and causal relations) and non-stationary data as well. GOLEM (Ng et al. 2020) lives at the interscction of score-based and gradient-based methods. It solves an unconstrained optimization problem where the objective function is given by a likelihood function (as in score-based methods), penalized by regularization terms for sparsity and acyclicity.
As alrcady mentioned, the concept of multiscale, non-stationary causal graphs is an understudied topic. Since none of the previous baseline models have been developed to infer causal graphs from data obeying to an underlying MN-DAG, our results provide information regarding the robustness of the previous algorithms with respect to the presence of multiple time scales and non-stationaritty In the following experiments, we use the code of MSCASTLE developed by D'Acunto et al. (2022): we exploit the implementations of DirectLiNGAM and GOLEM provided by gCastle" Zhang et al. 2021a whereas we resort to causalleart for the implementation of CD-NOD. The configuration for cach bascline is provided in Appendix F
Differently from MN-CASTLE,, baseline models are non-probabilistic. While our model provides an ap- proximate predictive posterior distribution over MN-DAGs, baseline models return a point estimate of an acyclic causal structure. In order to compare the algorithms, we retain all causal cocfficinnss identified by MN-CASTLE that are in modulus significantly greater than 0.1 at 99% level.
Performance in the estimation of the adjacency tensor. Retrieving the adjacency tensor means identifying the presence of causal relations disregarding their intensity. Appendix G provides example of the evolving causal relations inferred by MN-CASTLE, while Appendix H gives insights concerning the estimate of the non-stationarity parameter. Figure 7 refers to the first configuration described above, and shows the F1 score (the higher the better) for the considered models. The definition of the considered metries are given in Appendix I Appendix J also reports the values for additional evaluation scores and the fraction of undirected edges for MN-CASTLE.
Given a (4) setting, for each model we have 20 values of Fl. We use the box plot in order to visualize the inter-quartile range (IQR). In addition, within the box plot we overlay the F1 scores attained for each of
http:///gghhuuuccccc/msssss////s https://github.. phtp//mmppii///ausaalaaaaaa
14

---

Published in Transactions on Machine Learning Research (10/2023)
09
1.0-
0.5-
0.6-
0.4-
CASTLE
MECASTLE
0.3-
GOLEM
DINCLPNDA
0.
0.5
0.5
0.5
0.9
Figure 7: The figure depicts the performances of the considered methods in the retrieval of the adjacency tensor, according to F1 score. Higher F1 indicates better performance. Each model associated to a different color. For every (T, μi) setting and every model, we represent the values attained over the 20 synthetic datasets through box plot. Within the latter, we overlay the performances obtained for cach dataset (points).
the 20 datasets, plotted as points. For each value of т, we provide the performances of each algorithm as l varies. In case μ & O, we return the performance of GOLEM as well. In particular, because < is global, we replicate the causal structure retrieved by the latter for each time scale. This way, we obtain an additional baseline method also for multiscale datasets.
Overall, MN-CASTLE outperforms the bascline models in cach case. When μ 0, the performances of MSCASTLE and GOLEM is very similar to that of our method. In contrast, μ increases, we observe that the gap between MN-CASTLE and the other models widens and that, in general, the two multiscale models outperform GOLEM. In addition, MN-CASTLE performance remains stable as 7 increases, and when j increases the IQR tightens.
Figure 8a depicts the performances of the models as we vary the density of the underlying MN-DAG (second setting). Since here we use μ 0.5, we only retain MN-CASTLE, MSCASTLE and GOLEM. MN-CASTLE outperforms the baseline models along 8. In addition, the larger 8, the better the performances of our model.
Figure 8b provides the performances of the methods as we vary the size of the underlying MN-DAG (third setting). Since f = 0.5, we only compare MN-CASTLE, MSCASTLE and GOLEM also here. MN-CASTLE consistently outperforms the baseline models along N. The overall downtrend in the performances of the considered methods stems from the growth of the dimensionality of the problem while the number of obser- vations is kept fixed.
Last but not least, Figure &c depicts the performances of MN-CASTLE under model misspecification. The results prove that our model behaves good as the case with no violation in every scenario.
Performance in the estimation of 0. Figure provides results concerning the goodness of the inferred vector of scores i of PL distribution for the first experimental configuration, measured by normalized discounted cumulative gain (nDCG) at 3, w.r.t. the ground truth causal ordering x. Appendix I provides insights concerning the computation of this metric, whereas [ppendix J shows the results for Kendall-r statistics, Spearman's rank correlation and nDCG at 5. To represent the results, we use box plots. For each of the 20 synthetic datasets, generated according to specific pair (T, we obtain an estimated Ó. Then, we sample 103 causal orderings & from PT.()). Now, for each drawn causal ordering, we evaluate the monitored metric w.r.t. X. As vector of scores for bascline model, we use 0, where 0 ~ U(O, i i= I,...,N. N. Afterwards, we obtain 103 random causal orderings ス by sampling from the PL distribution parameterized by 0. As for MN-CASTLE, built by using 2 x 10a points. Overall, according to the monitored metric, MN-CASTLE outperforms the baseline model. In addition, the performances do not deteriorate as t grows and improve as increases.
evaluate the metric w.r.t. <. Therefore, for each model, every box plot is
We apply the same approach to evaluate the sensitivity of the estimation accuracy for 0 with respect to the density and number of nodes of the underlying MN-DAG and under generative model misspecification. Figure 10 shows the results obtained on the synthetic data generated according the second, third, and fourth experimental settings, described above.
15

---

Published in Transactions on Machine Learning Research (10/2023)
0.-
0.-
MSEASTTE
COLEM
0.2
05
0.75
(a) F1 6
0.5
0.-
0.2
N
(b) FI
1.0
0,5
0.4
vickrion
Uniomm
Subesmpiig
(c) under model misspeciifcattin
Figure 8: The figure depicts the performances of the considered methods in the retrieval of the adjacency tensor, according to F1 score, along (a)) 5 and (b) N, and (ch under generative model misspecification. In the latter setting, to ease the comparison also report the performances for the case without violation. Higher F1 indicates better performance. Each model is associated to different color. We represent the values attained over the 20 synthetic datasets through box plots, where we overlay the performances obtained for each dataset (points).
0.0
0.5
0.9
1.0
Model
02-
MN-CASTLE
Random
0.0
0.0
0.5
0.0 0.5
05
Figure The figure depicts box plots along with quartiles reference lines (dashed lines) for normalized discounted cumulative gain (nDCG) 3. MN-CASTLE is given in blue while random baseline model in orange. For every dataset generated according to given (r, ()) setting (i) sample 1 x 103 causal orderings 2 PL(O), where 0 is the estimated vector of scores; (ii) draw 1 x 10 random causal orderings evaluate nDCG@3 by using the sampled
R PT)()), where 0, (7O0N), i 1,..., N. Afterwards, causal orderings and x for both models. Thus, each box plot is made by 2 x 10 points.
16

---

Published in Transactions on Machine Learning Research (10/2023)
0.8-
nDC 61)
0.6-
0.2
MN-CASTLE
Random
0.0
0.25
0.5
0.75
10
Laplace
Uniform Violation
Subsampling
Figure 10: The figure depicts box plots for the normalized discounted cumulative gain (nDCG) at 0.6.V. MN-CASTLE is given in blue while random baseline model in orange. Box plots the left refer to the second experimental configuration, i.e, when we vary ó while keeping fixed the values of the others parameters, as described above. Box plots on the center concern the third experimental setting, where we study the sensitivity of the estimation accuracy w.r.t. the network size N. Box plots on the right relate the fourth experimental sctting, where we study MN-CASTLE under generative model misspecification. Notice that in the left and right plots N 5
Overall, the accuracy of MN-CASTLE in retrieving 0 grows along with 6: the IQR of the monitored metric reach higher values. In addition, the performance of MN-CASTLE in recovering the causal ordering is high and shows no dependence on N.. Notice that the nDCG score depends on the value of N. Thus, the comparison is meaningful because considers the same fraction of nodes for each combination. Finally, the performance MN-CASTLE does not degrade under model violations.
8 Analysis of Natural Gas Prices in the US Market
In this section, we examine the key drivers of natural gas prices in the US market during the period spanning from January 1, 2018, to December 31, 2022. Our analysis considers several variables, including the price of natural gas (NG), crude oil (CO), deviations in gas storage (SD), rig counts targeting gas (RC), deviations from scasonal average values of gas consumption for cooling (CDD) and heating (HDD) environments, the crack spread between heating oil and crude oil (CS), and the economic uncertainty index (UI, Baker et. al., 2016). We collected the data on weekly basis, and we analyze time scales ranging from 1 to 5, corresponding to resolutions ranging from 2-4 (scale 1) to 32-64 (scale 5) weeks. Please refer to Appendix K for further details on data sources and the pre-processing of the time series data.
The algorithms used for this dataset are GOLEM, MSCASTLE, and MN-CASTLE. GOLEM identifies stationary instantancous interactions between variables, meaning relationships that occur at frequency higher than weekly and remain constant over time. It does not establish any causal relationship with NG, as the only interaction detected is SD-HDD.
MSCASTLE is capable of detecting causal interactions between the time series on the 5 considered time scales. However, these interactions are assumed to be stationary over time. Figure I1 shows the multiscale causal networks learned by MSCASTLE. The edges associated with a positive causal coefficient. are repre- sented in green, while those with negative coefficient are represented in pink. The thickness of the edge is directly proportional to the absolute value of the causal coefficient. Unlike GOLEM, the multiscale analysis allows MSCASTLE to detect interactions that occur at longer time resolutions. The method suggests that NG causes CO on both scale 1 and 2, with positive causal cocfficient. While it is known from the literature that the price of CO drives the price of NG due to the substitution processes of NG with petroleum prod- ucts, the NG-CO relationship represents a novel element difficult to justify in light of the interchangeable relationships between CO and NG. Like GOLEM, MSCASTLE detects negative causal interaction from SD to HDD, which occurs on the first, third and fourth scales. On the contrary, on the second scale the relationship reversed. Overall, the interactions between SD and HDD are the larger in magnitude.
17

---

Published in Transactions on Machine Learning Research (10/2023)
C0
CO
U1
HDD
HDD
.RC
CDD
NG
CS
C$
Storage
Scale
Crude Oll:
Uncertainty Indea: Ul
HDD
HD0
Heating days
RC
Rig count: RC
COO'
Cooling days norm deviation: CDD
Crack spread: CS
Figure 11: The figure depicts the multiscale DAG retrieved by MSCASTLE
SD
U1
HDD
RC
CDD
Causal ordering x
2.5
Scale
0.0
Storage deviation:
-2.
Crude Oil: CO
2.5
Scale
0.0
Uncertainty Index: Ul
2.5
Heating days deviation:
2.3
Scale
0.0
Rig count: RC
2.5
2.5
Cooling days deviation: CDD
NG
0.0
2.5
Natural Gas: NG
Crack spread: CS
NG
0.0-
-35
Figure 12: The figure depicts the relationships inferred by MN-CASTLE, which involve NG at the considered scales, and the estimated causal ordering -
18

---

Published in Transactions on Machine Learning Research (10/2023)
MN-CASTLE can detect causal relationships between variables at different temporal resolutions, allowing for more nuanced insights into how the variables interact over time. The behavior of the causal coefficients learned by MN-CASTLE is illustrated in Figure 12 For readability, here we only report the coefficients involving NG, which is the focus of our analysis. The full MN-DAG is given in Appendix L Overall, the information provided by our method is richer and more detailed than the baselines. This is due to the ability of our method to track the temporal evolution of causal connections and detect relationships that are activated in a limited period of time and that may change sign. As a result, MN-CASTLE suggests denser causal structures, thus providing more comprehensive understanding of how the variables are related.
Interostingly, many of the causal relationships inferred by MN-CASTLE intensify during two specific periods: the early months of 2020 and 2022. These time windows correspond to significant events that put pressure on the energy market: the outbreak of COVID-19 and the Russian invasion of Ukraine. The fact that MN-CASTLE can detoct these changes causal relationships highlights its ability to capture the dynamic nature of the market and the impact of exogenous events.
Concerning the causal ordering, in Figure 12 we observe that SD, CO, and UI occupy the first three positions, whereas NG and CS occupy the last ones. Our method is the only one able to detect that the price of NG is causally influenced by scasonal factors, economic uncertainty, oil price, and deviations in gas reserves. These relationships occur at larger temporal resolutions of weeks, indicating that short-term fluctuations have a limited impact on NG prices. Our findings are consistent with the literature (Brown & Yucel 2008
Nick & Thoenes 2014 Ji et al. 2018), which highlights the relevance of economic and environmental factors in shaping energy markets. In line with the estimated causal ordering, MN-CASTLE does not propose the NG-CO interaction returned by MSCASTLE. Furthermore, the absence of interactions at lower scales is in line with the fact that GOLEM does not detect any connection regarding gas and underlines the superiority of multiscale methods. Regarding the impact of SD on NG, we see that at scales 2 and 3, an increase/ decrease in SD causes a decrease/ decrease/increase finceasee in NG prices, thus suggesting that supply shocks in the mid term can affect NG pricos. On the contrary, on scale 5, we observe a positive relationship around the years 2018-2019. This relationship can be explained by the fact that the demand for NG increased in that period (also due to exports) thus causing an increase in the price of NG, while gas reserves increased due to the record production of NG in the US.
MN-CASTLE detects the causal interaction of CO-NG at different time scales and highlights changes in sign over the analyzed period. Specifically, scales 3 and 4 show that in 2021, the increase in CO prices causes a corresponding increase in NG prices. However, in 2022, scales 2 and 4 indicate a negative relationship, possibly due to geopolitical tensions and speculative phenomena.
In addition, we observe that the prices of NG tend to be higher during periods with greater deviations HDD and CDD, which represent unusual weather events.
9 Conclusions and Future Research Directions
This paper deals with multiscale non-stationary causal analysis, filling a gap in the literature. Indeed, the bulk of previous work assumes that the only relevant temporal resolution for causal relations is the frequency of observed data. We drop such assumption. In addition, we also allow the causal relations to vary over time. Since general there no prior knowledge about the relevant time scales of causal interactions nor about their temporal dependencies, the proposed framework of MN-DAGs represents an important step in such direction.
Generative model. We propose probabilistic model to generate time series data obeying to underlying MN-DAG, in accordance with the specified values for multiscale and non-stationary features, l and т, respectively. Our model leverages the well established mathemattcca theory of multivariate locally stationary wavelet processes and linear structural equation model. The causal ordering is modeled by means of Plackett- Luce distribution while the causal interactions evolve over time according to the specified kernel of a Gaussian process. Statistical analysis of generated data proves the exposed model to be able to reproduce well-known features of time series. Therefore, it represents a suitable framework for testing the robustness of causal
19

---

Published in Transactions on Machine Learning Research (10/2023)
structure learning methodologies on datasets generated from different points of the (7, ()-qadaatt shown in Figure 1
We stress the importance of providing both researchers and practitioners with synthetic data generators capable of replicating phenomena characterizing data from different application domains.
Future work should aim to overcome some limitations related to the framework adopted to manage different time resolutions and the modeling of the causal tensor. In particular, multivariate locally stationary wavelet processes formulation relies upon wavelets, that are known to suffer from limited joint time-frequency reso- lution (Heisenberg uncertainty principle). Indeed, wavelets divide the frequency space into non-overlapping bands, Le., octave bands. Furthermore, since the auto/cross-con /cross-coorreation structure of generated data depends on both the power spectrum decomposition across temporal scales and the auto-correlation wavelet, the usage of diverse wavelet families might lead to different results. Then, the usage of alternative methods to wavelet transform might improve the proposed generative model. With regards the causal tensor, structural breaks such as sudden leletion/addition addition of causal edges might be added within [Equatio (1)
From theoretical point of view, an interesting research direction is to study the assumptions that make the model described in Eq. (2) idcntifiable. Even though some class of lincar structural equation models have been proved identifiable under different types of restrictions (Shimizu et al. 2006 Peters & Bühhlannn, 2014 Loh & Bühlmann 2014: Park & Kim 2020), the case of MN-DAG needs to be carefully investigated.
Indeed, the presence of the non-decimated wavelet transform: the unobservability of the contributions to the
process coming from each time resolution; the linearity of the model in the frequency domain are some of the points that distinguish the MN-DAG case from those currently studied.
Bayesian causal structure learning method. In addition, we expose a Bayesian method for learning MN-DAGs from time series data, termed MN-CASTLE, The latter relies upon observed time series data and an estimate for the inverse power spectrum at each scale level. We implement the latter by using a two-step approach. In the first step we optimize w.r.t. the Plackett-Luce vector of scores 9, by using the values of time series at time f. Then, we keep the causal ordering fixed to the mode of the Plackett-Luuce distribution, i.e., <0, and we estimate the rest of variational parameters related to the cansal coefficient tensor by exploiting the provided estimation for the inverse power spectrum.
Our findings show that MN-CASTLE compares favorably to baseline models in the retrieval of the adjacency tensor of the causal graph. We test the models synthetic datasets generated according to different (T, p)) configurations, from the single-scale stationary to highly multiscale non-stationary case. We observe that the performance of MN-CASTLE, depicted in Figure 7| is not sensitive to the value of non-stationarity parameter. On the contrary, the growth of the multiscale parameter is associated with an improvement in the quality of the results returned by our method since the IQR tightens. The improvement in performances along j also shown in Figure 9, that concerns the goodness of the estimated vector of scores for Plackett- Luce distribution. On one hand, we think that when non-stationarity and multiscale parameters are different from zero, MN-CASTLE might benefit of greater differences among time scries distributions. On the other hand, we believe that the large variance shown (especially in the single scale case) is effect due to the low cardinality of the edge set. In fact, even though the monitored metric is normalized, on average in the single scale case we only have five causal links. So, a single error weighs more. We emphasize that MN-CASTLE, being a fully Bayesian approach, by definition takes into account uncertainty. Consequently, we might sample MN-DAGs from the approximate posterior distribution in accordance with the confidence of the model.
Furthermore, study the behaviour of our model w.r.t. the density ö of the underlying MN-DAG. We observe that the performance of MN-CASTLE improves as ò increases and that our method outperforms the other models in all cases. This improvement is also manifest in the value of the metric used to evaluate the estimated causal ordering. We also provide supplementary results on additional synthetic data to lest the capabilities of the monitored methods when the MN-DAG size N increases, keeping the other parameters fixed. Also in this settings, MN-CASTLE outperforms the basclines. In addition, the ability of MN-CASTLE in estimating causal ordering does not deteriorate as N increases. Last but not least, MN-CASTLE keeps performing well under model misspecification, thus highlighting the robustness of our method.
20

---

Published in Transactions on Machine Learning Research (10/2023)
In case study, we have applied MN-CASTLE to analyze the drivers of natural gas in the US market and have compared the results of our method with those of GOLEM and MSCASTLE. While cannot determine the ground-truth, we have found that MN-CASTLE provides richer information than the baselines, as it can track the evolution of causal relationships across different scales over time. Our study has revealed that causal relationships have strengthened during the outbreak of COVTD-19 and at. the beginning of the Russian invasion of Ukraine. Additionally, by accounting for non-stationarity, MN-CASTLE has detected more relationships than the baselines. Furthermore, our method is the only one capable of detecting the causal impact of seasonal factors, economic uncertainty, oil prices, and deviations in gas storage on natural gas prices, which are crucial drivers in the Economics literature.
References
Yashas Annadani, Jonas Rothfuss, Alexandre Lacoste, Nino Scherrer, Anirudh Goyal, Yoshua Bengio, and Stefan Bauer. Variational causal networks: Approximate Bayesian inference over causal structures. arviv preprint Xrx::22066. 0e21106006633 07635, 2021.
Scott R Baker, Nicholas Bloom, and Steven J Davis. Mcasuring cconomic policy uncertainty. The Quarterly Journal of Economics, 131 (4):1593 1636, 2016.
Matthias Bauer, Mark van der Wilk, and Carl Edward Rasmussen. Understanding probabilistic sparse Gaussian process approximations. Advances in Neural Information Processing Systems, 29, 2016.
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 566, 2010.
Eli Bingham, Jonathan P Chen, Martin Jankowiak, Fritz Obermeyer, Neeraj Pradhan, Theofanis Karaletsos, Rohit Singh, Paul Szerlip, Paul Horsfall, and Noah D Goodman. Pyro: Deep universal probabilistic programming. Journal of Machine Learning Research, 90(1):973-978 2019.
Christopher M Bishop and Nasser M Nasrabadi. Pottern recognition and machine learning, volume 4(4).
Springer, New York, New York, 2006.
David M Blei, Alp Kucukelbir, and Jon D McAuliffe. Variational inference: A review for statisticians. Journal of the Americon Statistical Association, 112(518):859 877, 2017.
Stephen PA Brown and Mine K Yucel. What drives natural gas prices? The Energy Journal, 29(2), 2008.
Peter Bühhmann, Jonas Peters, Jan Ernest, et al. CAM: Causal additive models, high-dimensional order search and penalized regression. The Arnals of Statistics, 42(6):2526 2556, 2014.
Bertrand Charpentier, Simon Kibler, and Stephan Günnemann. Differentiable DAG sampling. arXiv preprinl arXi::2203.0850 2022.
David Maxwell Chickering. Optimal structure identification with greedy search. Journal of Machine Learning Research, 3 5504 2002.
Chris Cundy, Aditya Grover, and Stefano Ermon. Bed nets: Scalable variational approaches for Bayesian causal discovery. Advances in Neural Information Processing Systems, 34:7095 7110, 2021.
Gabriele D' 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. 8, 2021.
Gabricle D' D'uunto, Paolo Di Lorenzo, and Sergio Barbarossa. Multiscale causal structure lcarning. CoRR, abs/2207.07900 079008, 2022. doi: 11.48550/arXiv22000 07908
URL https:://dci. 4855000r05i0.0 ar20
Rainer Dahlhaus. Fitting time series models to nonstationary processes. 1-37, 1997.
The Annals of Statistics, 25(1):
21

---

Published in Transactions on Machine Learning Research (10/2023)
Rainer Dahlhaus. Graphical interaction models for multivariate time series. Metrika, 51:157-172, 2000.
David A Dickey and Wayne A Fuller. Distribution of the estimators for autorogressive time scries with a unit root. Journal of the American Slatistical Associalion, 74(366a):427 431, 1979.
Zhuanxin Ding and Clive WJ Granger. Modeling volatility persistence of speculative returns: A new ap- proach. Journal of Econometries, 73(1):185-22) 1996.
Idris A Eckley and Guy Nason. Efficient computation of the discrete autocorrelation wavelet inner product matrix. Statistics and Computing, 15(2):83-99) -2,, 2005.
Piotr Fryzlewicz, Sébastien Van Bellegem, and Rainer Von Sachs. Forecasting non-stationary time series by wavelet. process modelling. Annals of the Institute of Slatistical Mathematics, 55(4):737-70) 2003.
Artyom Gadetsky, Kirill Struminsky, Christopher Robinson, Novi Quadrianto, and Dmnity Vetrov. Low- variance black-box gradient estimates for the Plackett-Tuce distribution. In Proceedings of the AAAT Conference on Artificial Intelligence, volume 34(06), pp. 10126 10135, 2020.
Jacob Gardner, Geoff Pleiss, Kilian Q Weinberger, David Bindel, and Andrew G Wilson. Gpytorch: Black- box matrix-matrix Gaussian process inference with GPU acceleration. Advances in Newral Information Processing Systems, 31, 2018.
AmirEmad Ghassami, Negar Kiyavash, Biwei Huang, and Kun Zhang. Multi-domain causal structure learn- ing in linear systems. Advances in Neural Information Processing Systems, 31, 2018.
Clark Glymour, Kun Zhang, and Peter Spirtes. Review of causal discovery methods based on graphical models. Frontiers in Genetics, 10:524, 2019.
Mingming Gong, Kun Zhang, Bernhard Schölkopf, Dacheng Tao, and Philipp Geiger. Discovering temporal causal relations from subsampled data. In International Conference on Mochine Learning, pp. 1898-18888 1906. PMLR, 2015.
Tan Goodfellow, Yoshua Bengio, and Aaron Courville. Deep Learning. MIT Press, Cambridge, MA, USA, 2016. http://wwr http://www.deepleaaae www.oee. org
David Heckerman, Dan Geiger, and David M Chickering. Learning Bayesian networks: The combination of knowledge and statistical data. Machine learning, 20(3):197 243, 1995.
James Hensman, Alexander Matthews, and Zoubin Ghahramani. Scalable variational Gaussian process classification. In Artificial Intelligence and Statistics, pp. 351-360. PMLR, 2015.
Matthew D Hoffman, David M Blei, Chong Wang, and John Paisley. Stochastic variational inference. Journal of Machine Learning Research, 2013.
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, 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 Ramsey, Ruben Sanchez-Romero, Clark Glymour, and Bern- hard Schölkopf. Causal discovery from neterogeneouss/oonstttiion data. Journal of Machine Learning Research, 21(89)11-53, 2020.
Carlos M Jarque and Anil K Bera. A test for normality of observations and regression residuals. International Statistical Review/Rever Internationale de Statistique, pp. 163-172, 1987.
Qiang Ji, Hai- Ying Zhang, and Jiang- Bo Geng. What drives natural gas prices in the United States? A directed acyclic graph approach. Energy Economics, 69:79 88, 2018.
22

---

Published in Transactions on Machine Learning Research (10/2023)
Nan Rosemary Ke, Olexa Bilaniuk, Anirudh Goyal, Stefan Bauer, Hugo Larochelle, Bernhard Schölkopf, Michael C Mozer, Chris Pal, and Yoshua Bengio. Learning neural causal models from unknown interven- tions. arxiv preprint orXiv1910.00077. 2019.
Diederik P Kingma and Jimmy Ba. 1.X.v:112.66988 2014.
Adam: A method for stochastic optimization. arXiv preprint
Diederik P Kingma and Max Welling. Auto-Encoding Variational Bayes. arXiv preprint arXiw:1312.6114, 2013.
Po-Ling Loh and Peter Bühlmann. High-dinensional learning of linear causal networks via inverse covariance estimation. Journal of Machine Learning Research, 15(1) ):3065-331000 2014.
Lars Lorch, Scott Sussex, Jonas Rothfuss, Andreas Krause, and Bernhard Schölkopf. Amortized inference for causal structure learning. orXiv preprint arXiv:2205.122 12934, 2022.
R. Duncan Luce. Individual choice behavior: A theoretical analysis. John Wiley & Sons, New York, 1959.
Benoit Mandelbrot. The variation of some other speculative prices. The Journal of Business, 40(4):393 413, 1967.
Andriy Mnih and Karol Gregor. Neural variational inference and learning in belief networks. In International Conference on Machine Learning, pp. 1791 1799. PMLR, 2014.
Guy P Nason and Bernard W Silverman. The stationary wavelet transform and some statistical applications. In Wavelels and statistics, pp. 281 299. Springer, 1995.
Guy P Nason, Rainer Von Sachs, and Gerald Kroisandt. Wavelet processes and adaptive estimation of the evolutionary wavelet spectrum. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 62(2))271-292 2000.
Tgnavier 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.
Ignavier Ng, Shengyu Zhu, Zhuangyan Fang, Haoyang Li, Zhitang Chen, and Jun Wang. Masked gradient- based causal structure learning. In Proceedings of the 2022 SIAM International Conference on Data Mining (SDM), pp. 12--432. SIAM, 2022.
Sebastian Nick and Stefan Thoenes. What drives natural gas prices? A structural VAR approach. Energy Economics, 45:517-527, 2014.
Gunwoong Park and Youngwhan Kim. Tdentifiability of Gaussian linear structural equation models with homogeneous and heterogeneous error variances. Journal of the Korean Statistical Society, 49(1):276 292, 2020.
Timothy Park, Idris A Eckley, and Hernando C Ombao. Estimating time-evolving partial coherence between signals via multivariate locally stationary wavelet processes. IEEE Transactions on Signal Processing, 62 (20):5240-5200, 2014.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zemiin Lin, Natalia Gimelshein, Luca Antiga, et al. Pytorch: An imperative style, high-performance deep learning library. Advances in Neural Information Processing Systems, 32, 2019.
Judea Pearl. Causality. Cambridge University Press, Cambridge, UK, 2009.
Ronan Perry, Julius von Kügelgen, and Bernhard Schölkopf. Causal discovery in heterogencous environments under the sparse mechanism shift hypothesis. arXiv preprint arXiw:2206.0001: 020013 2022.
Jonas Peters and Peter Büihmann. Identifiability of Gaussian structural equation models with equal error variances. Biometrika, 101 (01(1)211 (1::):222222 228, 2014.
23

---

Published in Transactions on Machine Learning Research (10/2023)
Jonas Peters, Joris M Mooij, Dominik Janzing, and Bernhard Schölkopf. Causal discovery with continuous additive noise models. Journal of Machine Learning Research, 2014.
Robin L Plackett. The analysis of permutations. Journal of the Royal Statistical Society: Series (Applied Statistics), 24(2):193-20) 202, 1975.
Bechir Raggad. Time varying causal relationship between renewable energy consumption, oil prices and economic activity: New evidence from the United States. Resources Policy, 74:102422, 2021.
Jakob Runge, Sebastian Bathiany, Frik Bollt, Gustau Camps-Valls, Dim Coumou, Ethan Deyle, Clark Glymour, Marlene Kretschmer, Miguel D Mahecha, Jordi Muñoz-Mari, et al. Inferring causation from time series in Earth system sciences. Nature Communiicttioss, 10(1):1-13, 2019.
Robert B Schnabel and Elizabeth Eskow. A revised modified Cholesky factorization algorithm. SIAM Journal on Oplimizalion, ((4):1135-1144 1999.
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): 612-634, 2021.
Shohei Shimizu, Patrik O Hoyer, Aapo Hyvärinen, Antti Kerminen, and Michael Jordan. A linear non- Gaussian acyclic model for causal discovery. Journal of Machine Learning Research, 7(10), 2006.
Shohei Shimizu, Takanori Inazumi, Yasuhiro Sogawa, Aapo Hyvärinen, Yoshinobu Kawahara, Takashi Washio, Patrik O Hoyer, and Kenneth Bollen. DirecLLiNGAM: A direct method for learning a linear non-Gaussian structural equation model. Journal of Machine Learning Research, 12:1225 1-225-1488, 2011.
Le Song, Mladen Kolar, and Eric Xing. varying dynamic Bayesian networks. Advances in Neural Information Processing Systems, 22, 2009.
James C Spall. Introduction to stochastic search and optimization: Estimation, simulation, and control, volume 65. John Wiley & Sons, Hoboken, NJ, USA, 2005.
Peter Spirtes, Clark N Glymour, Richard Scheines, and David Heckerman. Causation, prediction, and search. MIT Press, Cambridge, MA, USA, 2000.
Eric V Strobl. Improved causal discovery from longitudinal data using a mixture of dags. In The 2019 ACM SIGKDD Workshop on Causal Discovery, pp. 100-133. PMLR, 2019.
Simon AC Taylor, Timothy Park, and Idris A Eckley. Multivariate locally stationary wavelet analysis with the mvLSW R. package. Journal of Statistical Software, 90:1-19, 2019.
Michalis K Titsias. Variational model selection for sparse Gaussian process regression. Report, University of Manchester, UK, 2009.
Ronald J Williams. Simple statistical
algorithms for connectionist reinforcement learning.
Machine Learning, 8(3):229-22 256, 1992.
Keli Zhang, Shengyu Zhu, Marcus Kalander, Tgnavier Ng, Junjian Ye, Zhitang Chen, and Lujia Pan. gCastle: A Python toolbox for causal discovery. arXiv preprint orXie:2111.15155 2021a.
Kun Zhang and Aapo Hyvärinen. Distinguishing causes from effects using nonlinear acyclic causal models. In Causality: Objectives and Assessment, pp. 157 164. PMLR, 2010.
Xun Zhang, Fengbin Lu, Rui Tao, and Shouyang Wang. The time-varying causal relationship between the Bitcoin market and internet attention. Financial Innovation, 7:1-19, 2021b.
24

---

Published in Transactions on Machine Learning Research (10/2023)
Appendix A Linear Structural Equation Models and DAGs
Mathematically, a DAG is formulated as SEM. Given a datasct X == (I1,...,, IN) of N random variables, SEM is collection of N structural assignments
Ti = Si(Pi, 2i), i == 1......, N.
where Pi represents the set of direct causes (parents) of node 24, is noise variable satisfying 2; 1 2j if * i, and fi(-) is generic functional form. In this paper, we focus on linear functional forms, therefore, by exploiting matrix form, the equation above becomes:
=
where C CeD N is the matrix of causal coefficients satisfying (i) ca =0 Vi (ii Gj #0 =
2, E Pú; (iii) diag(C") = O, Vn E N (acyclicity property). Since I - C is an invertible matrix (see Lemma Appendix F. we can rewrite the latter equation as
= MZ.
(4)
with M = -I- being mixing matrix. According to Equation (4) observed data is a mixing of independent latent noises. Here, causal relations are stationary, instantaneous and are supposed to occur at the frequency of observed data.
Appendix B Locally Stationary Wavelet Process
In the univariate case, locally stationary wavelet process (LSW, Nason et al. 2000 is a suitable modeling framework to represent a non-stationary process Xy of length T = 2 E N, by means of a triangular multiscale representation
+00
2T
Σ
Σ
by o k..
(5)
j=1
The building blocks of Equation (5)
are: (i) the random amplitude u/[4/T]2j.. composed by a time-varying
amplitude v, [k/]]] and normal noise variable 2j,4 such that cou(2 ) , where رر.ر' represents the Kronecker delta; (ii) discrete, real valued and compactly-supportee oscillatory functions t namely non-decimaled wavelels. At each time only some values contribute to xr [4], and the lime-dependence is managed by the index k. Local stationarity means that the statistical properties of the process vary slowly over time. This feature is essential order make learning possible (Nason et al. 2000). Within the LSW franework, local stationarity is formalized by means of a smoothness assumption concerning the time-varying amplitudes My [k:/T] (Fryzlewicz et al. 2003). Indeed, the latter quantity provides a measure of the time- dependent contribution to the variance
certain time scale level J, namely the evolutionary wavelet spectrum (EWS), defined as Sy[n] = |u,[,]] with v = k/T being the rescaled time (Dahlhaus) 1997). For a stationary process, FWS is constant vj < J. As an example, consider the MA(1) = T/V2(c[ c 1)). We obtain it by setting in Eqquatio (5) the following values for the previous components: (i) 2j,b cit] (fi) S, = if j = and zero otherwise; d = [1/V2, 1122.11/33 1/V2] as the Haar wavelet. Because S, constant and different from zero only for j = we obtain stationary amplitude wy [v] only for the first scale level. Then, it follows that
IT|
By tylk/T]zzAdgltt k]
1 1.6,6
1
V2 (e(t) - -elt - 1])
25

---

Published in Transactions on Machine Learning Research (10/2023)
Appendix C Multivariate Locally Stationary Wavelet Process
The MLSW framework generalizes locally stationary wovelet process (LSW, Nason et al. 2000 see Appendix
B to model N zero-mean processes Xr!4]
[|] =
each of length T,, as follows:
J Xrl4| Σ Σ Vylk/T/aa kj.
(6)
j=
In Equation (6) (1) [by[t=)] a set of non-decimated wavelets; (ii) {2,,}} a set of random vectors 2j.k (O, [xx);; (iii V, [h/T] E NNN is the transfer function matrix, assumed to be lower triangular and with entries being Lipschitz continuous functions associated with Lipschitz constants L, 1,f
[V) e {I,...,N), such that E, -. (.)m) < co. Local stationarity means that the statistical properties of the process vary slowly over time. This feature is essential in order to make learning possible (Nason et al., 2000). and within MLSW coincides with the Lipschitzianily assumption above. Here, the transfer function malrix V,[e), with v = k/T being the rescaled time (Dahhhaus 1997), provides measure of the local variance and cross-covariance between the processes at a certain time y and scale j, i.e. the local wovelet spectral matric (LWSM) S, [/] = V)[vVVj[r).
By construction, LWSM is symmetric and positive cach time and scale j. Within LWSM, diagonal elements S%,, [r] represents the spectra of of the processes, whereas S [v] provides the cross-spectra wee them. In addition, the local auto and cros-covariance functions, namely (v,(v,)) and Cam (v, () (with I being certain lag), admit formulation in terms of the LWSM (see Park et al. 2014 for further details).
Appendix D Overview of SVI
Stochastic variational inference (Hoffman et al. 2013 Kingma & Welling 2013 is algorithm that combines variational inference (VI, Blei et al. 2017]] and stochastic optimization (Spall, 2005). SVI approximates the posterior distribution of complex probabilistic models that involves hidden variables, and can handle large datasots. Consider a dataset x= [T of T ii.. samples of either a continuous or discrete variable x. Suppose that X generated according to latent coutinuous random variable z, The latter is governed by vector of parameters 8 endowed with a prior distribution p(8") i.e., z(6) ~ (z). Thus, we have data are generated according to conditional distribution, ie.,, x(i) Pa- (x z). Both the prior Pæ- (z) and the conditional distribution Pa- z) belong to parametric families of distributions ps(z) and ps(x z) whose PDFs are differentiable w.r.t. 9 and z. Our goal is to compute the likelihood of the hidden variable given the observations, i..., the posterior
(x()
ps(z x)
(7)
J (x, z) dz
Since the denominator of Equation (7) also known evidence, is usually intractable to compute, well- known solution is to approximate the target posterior. Within approximate posterior inference methodolo- gies, VT casts learning as an optimization problem. More in details, VI involves the introduction of family [z] x), parameterized by a variational parameters 6. Then, VI optimzes those parameters to find 90- (z x), i.c., the member of the variational distributions family that is closest to the posterior distribution. Here closeness is measured according to Kullback-Leibler divergence (KL).
of variational distributions 4
The objective of SVI the evidence lower bound (ELBO), that is equal to the negative KL divergence up to term that does not depend on
ELBO- Pooleex)looppxxx, 2) log (z x)]
DDDuzqgslz (1) (z x ())) + log (x)
DKI(go(Z x (0) [z)) + (x)))) [log (x [z)].
Since KL is nen-negative measure of closeness between distributions, then log 2 ELBO for all B and sb Therefore, the maximization of the F.L.OO is equivalent to the minimization of the distance between
26

---

Published in Transactions on Machine Learning Research (10/2023)
(z( and (x z). Observations x0 are conditionally independent given the latent, thus the log likelihood
term in Equation (8) can be written as
Σ log p(
2) R
T T
Σ
logp(x 22),
iEZT
where Ir is a set of indexes size T' 2 T. One way to subsample indexes is, for example, to randomly select T" data points among the obscrvations Thus, in case large datasots, we can run SVI while exploiting mini-batch oplimization.
In order to compute the gradient of the ELBO w.r.t. Φ SVI relies upon the reparameterization trick. The continuous random variable z can be expressed in terms of deterministic function z = _4,c, where ~ q(E) is independent of z. This procedure is useful to move all the dependence on & inside the expectation operator
ta (() = (ogess x
where f0(z) represents a general cost function. Now, the gradient can computed as
V V Voleatolfoloo(:x(0)) Egte) V Vofoloo(ex
1
f(os(cl.))x
where L is the number of samples per data point. Then, obtain unbiased estimate of the gradient by means of Monte-Carlo estimates of this expectation.
Appendix E Proofs
Lemma Appendix E.1. The inverse of T - C erists and consists of a finile sum of powers of C.
Proof. To prove invertibility of C- E NxN let us rewrite C = P'ĆP. Here, E NxN is a permutation matrix entailed by the causal ordering <, such that Pw'n = 1 ifT the node X,, occurs at position n' within x, and Č is strictly lower triangular matrix, computed by ordering the rows of C according Now, for permutation matrices it holds P-1 P'. In addition, since strictly lower triangular matrices are nilpotent, there exists an integer N such that C 0, Vn N. Then it follows that C is similar to C and, consequently, nilpotent too:
CN (P'ČP)
(P- (P-'ČP)
(P- CPP (P-(PP
P ")(PP=))((PP (PP (PP-))
p
= 0.
At this point, exploiting the geometric series representation (nilpotent matrices have eigenvalues equal to zero and then are convergent), we have that
(I (1- C)
n==
N
m=0
27

---

Published in Transactions on Machine Learning Research (10/2023)
Therefore, the inverse exists and is given by finite sum of powers of C.
Lemma 5.1. The transfer function matrix is permuled lower Iriangular malria, M2 = P'(T- P'( where P E RN*N is permatation matria such that Pu'n 1 iff the node Xn occurs at position n' within the causal ordering 3, and C,(0] is strictly lower triangulor matria of causal coefficients.
Proof. Starting from the representation of C,(4] = [uj, we have:
M2 (1 [PP)
(P- 1PP)
(P'P - P'C,((P))-
(P'(I- - C,(p1)P))
= -P- - [1,1-
where (I- (I- admits a representation in terms of the geometric series for Lemma Appendix E.1 which in this case consists in a sum of lower triangular matrices.
Lemma 5.2. The local wavelet spectral matriz and its inverse are given by S,(al = P'[I P'(I-Č,lvl) C,lu))-'P and Oyln] = P'(I -
Proof. For vellev invertible matrix A the Gramian AA' is semi-positive definite, hence Sylr] M/[jMM,lu] analogously to Definition 2 in Park al. (2014).
Then, by pluggin in the expression of My[t] given by Lemina 5.1
get:
s = (P'( - -1
P' Č, Č,-)-
where we exploit the fact properties P' = P-1 By following the same rationale and remembering that (BA)-1=A* --B-11 we obtain O,[r] = P - -
Proposition 5.3. The spectral properties of the process XT are independent of both the causal ordering and the order in which the process dimensions are observed.
Proof. Let us consider a causal matrix C,(y] and its causally ordered form Č,(z], such that Cj[y] P'Č Now, let us consider Sy[u] = (T C,(,)) - and Ö,(r] (T C,(+))(( - C2(1)) From Lemma 5.2 it directly follows that S (I- CIM))- C,(+))-*
P'S and Oylv]
(I-C, (I)))()))( (I = ['j, This proves independence from the causal ordering.
Now, let consider a process X, admitting representation of the form of Equation (2) Consider also a permuted version Xr, i.e, _e,,X== = QX with Q being an arbitrary permutation matrix of size x N. By multiplying both sides of Equation (2) by Q, we get that M, QM, (p]. Hence, we have that
S [Sy[r]Q' and [, = Q0,[r]Q', where Sy[²] = My[MM [p]' and O, [r] = S,Pp]-1 This proves
independence from the ordering of the dimensions of the process.
Appendix F Models Configuration
Below we report the models hyper-parameters used during the test phase:
28

---

Published in Transactions on Machine Learning Research (10/2023)
MN-CASTLE: fraction of inducing points equal to 64%; K KRBE: in case (the estimated Ō, is constant)) we use prior for Ax normal N(1. 10 1. 10-3); = 0.05; number of iterations iter= 6. x 102 with 10 particles;
MSCASTLE: & penalty parameter À 1. 10-1; pruning threshold 7 = 5. x 10 Daubechies wavelet with filter length equal to 2; maximum value for dagness function Ato = X 10-8;
GOLEM: pruning threshold 7 5. 10-2; number of iterations iter= 1. 10";
DirectLiNGAM: pruning threshold 7 5. x 10-2;
CDNOD: independence test = Fisher's 7; significance level 95%.
Appendix G Evolution over Time of the Estimated Causal Relations
We apply MN-CASTLE, over synthetic dataset constituted by N = time series of length T = 100 each. To generate the data, we use the exposed probabilistic model over MN-DAGs, where we set = å = 0.5 and we use as K = KRAP. In this case, we obtain = 3 scale levels. The configuration of MN-CASTLE is the same as Appendix F
Scale
Scale
Scale
0.0
2.5
0.0
0.0
-2.5
0.0
2.5
Estimated Actual
-25
0
10
17
Figure 13: The figure depicts the evolution over time (v= L/T) of estimated causal coefficients (blue) vs the ground truth latent coefficients (dashed orange), for the three temporal resolutions. Shaded bands refer to 99% confidence level.
F Figure 13 depicts the estimated causal relations and their evolution over time. LIN-CASLLE correctly tracks the behaviour of latent causal coefficients in all cases. As given in [Section 6 in the second inference step we use the mode of PL distribution 2u to mask the distribution of hidden functions Consequently, only those relations that conform to the estimated causal ordering show tight 99% Cls.
Appendix H Estimate of т
the results below we obtain the estimate by means of the estimated GP kernel lengthscale, i.c. /ARBF. In our experimental assessment, data are generated by using a RBF kernel. Then we can evaluate the goodness of A. However, this assumption might be too restrictive real-world scenarios, and we might want to use combination of kernels in the inference procedure. Therefore, the capability of computing lau might be impaired since might not have single lengthscale resulting from the kernel combination even we use kernel that does not involve any lengthscale in its formulation. Having said that, below we provide some insights concerning th estimation of the non-stationarity parameter т, coming from all the four experimental settings (see Section 7.2
Figure 14 depiets the inferred values for the non-stationarity parameter the first experimental setting. Here red dashed lines refer to the ground truth value of 7. Our model tends slightly overestimate the non-stationarity parameter. As the value of p increases, the estimate approaches the ground truth.
29

---

Published in Transactions on Machine Learning Research (10/2023)
59
33-
50
022
0.0
0.5
0.2
0.0
0.5
2b
Figure 14: The figure illustrates the box-plots concerning the estimated values e for the non-stationarity parameter. Red dashed lines refer to the ground truth value of T.
0.5
1.0
0.8
0.4
11.2
Laplace
Uniform Subsampling
Violation
Figure 15: The figure illustrates the box plots concerning the estimated values for the non-stationarity parameter. Red dashed lines refer to the ground truth value т. On the left, we have the results for the second experimental setting, on the center are those for the third, and on the right are those for the fourth.
In addition, Figure 15
provides the results for the second, third and fourth settings. Overall, we do not ap- preciate any trend along 6, N, and generalive model violations. Indeed, MN-CASTLE tends to overestimate the non-stationarity as in the first setting.
Appendix Definitions of the Performance Metrics
In this section we describe the metrics used to evaluate the goodness of the estimated adjacency tensor of the causal graph and the retrieved cansal ordering.
Adjacency. For the predicted adjacency tensor, monitor both accuracy and structural scores.
With regards to accuracy measures, we look at the true positive rate (TPR, recall), the false discovery rate (FDR, 1-precision) and the Fl-score. The first defined as TP/P, where TP the number predicted odges that exist in the ground truth with the same direction and P (condition positive) the number of links in the ground truth. The second given by FP/(FP+ TP). Here, FP is the number of edges that do not exist in the skeleton of the ground truth, i.e. in the undirected adjaceney. Finally, the F1-score is computed as the harmonic mean between TPR and FDR (precision).
Concerning structural metrics, first we consider the structural Hamming distance (SHD), that represents the number modifications (added, removed, reversed edges) needed to retrieve the ground truth starting from the estimated network. Then, we also monitor the ratio between the number of predicted edges and the condition positive, given by NNZ/P, where NNZ represents the sum of dirccted (D) and undirected (U) estimated edges. Finally, we have the fraction of predicted undirected edges, computed as FU U/NNZ.
Causal Ordering. To compare the estimated causal ordering with the ground truth, we consider three metrics able to provide measure of the association strength between two rankings.
30

---

Published in Transactions on Machine Learning Research (10/2023)
First, we look at Kendall-T, which is mcasure of ordinal correspondence between two rankings, bounded between -1 (low correspondence) and 1 (strong correspondence). Given two orderings 2 and x, the statistics is defined as:
Kendall-r = (P- (P-Q)
where here P is the number of concordant pairs, Q the number of discordant pairs, T the number of ties only in 2, and U the number of ties only in <;
Second, we employ & measure of ranking quality widely applied in information retrieval, the normalized discounted cumulolive goin (nDCG). Consider a ground truth ordering < of length N and suppose to associate items with descending scores s, from N to 1. Then, consider an other ordering 2 over the same set of elements in x. Now, define the discounted cumulative gain (DCG) as:
DCG
loga + 1)
and let the ideal discounted currulative gain (IDCG) to be the DCG of <. Therefore, the nDCG is defined as the ratio by the DCG and the IDCG. This score is bounded between the nDCG of the worst ordering of scores s, i.e., sorted in ascending order, and 1. In our analysis we use a min-max scaling map nDCG the unit interval. To evaluate the capability of a method in providing high-score items at first positions k. we compute the nDCG@k by considering only the first l clements of R.
Finally, we consider Spearman's rank correlation Ps, that provides a non non-parametric correlation coefficient between two series. Here, differently from Pearson correlation, data is not assumed to be normally dis- tributed. Thus, as Kendall-r, this statistics is bounded between -1 and 1. Since in our case we have two score vectors, namely 2 and K, made by distinct values, this metrics can be computed as:
PS 1 6 62 N(N2 - 1)
Appendix J Additional Monitored Metrics
In this section we provide additional analysis to better understand the behaviour of the considered methods when (i) we navigate the (1)quddat keeping the other parameters fixed, (ii) vary the density of the MN-DAG at point in the quadrant, (iii) change the size of the MN-DAG at point in the quadrant, and (iv) we violate the assumptions of the generative model. Figure 16a where SHD has been normalized by the number of edges present in the ground truth, refers to the first experimental setting and shows that MN-CASTLE provides the better performances in all scttings. Additionally, Figure 16b shows that MN-CASTLE reduces the number of false discoveries returned by baseline models, especially when μu / O. On the contrary, the single-scale stationary case, best FDR values are provided by GOLEM, Figure 16c tells that MN-CASTLE the best in the retrieval of true positives in almost all cases. Furthermore, the estimated to true edge set ratio given in Figure 16d shows that our model is the more accurate in the number of causal connections. Finally, we plot the fractions of undirected connections given by our model to check aciclicity of the retrieved structure.
With regards to causal ordering estimation, Figure 17 provides results also for Kendall-r ( Figuue 17a) Spearman's 's rank (Figure 17b) and nDCG@5 ( Figure 17ch statistics. Here, the methodology used the same as in Section 7.2 According these metrics, MN-CASTLE outperforms the random bascline model in all cases. In addition, the values of the monitored metrics do not lower as grows and improve as μ increases.
Hence, Figures 18 to 20 depict the resulting values for the same metrics above, obtained in the second, third, and fourth experimental settings, respectively. Overall, when the density of the network grows, the metrics improve for all methods. Also, MN-CASTLE provides the best performance for all values of 8. When we vary the network size, our method outperforms the baselines as well. The overall worsening in method's performance is due to the fact that here the dimensionality of the problem grows, while the number of observations is kept fixed, T 100. In addition, the performance of our method remains good even under generative model misspecification.
31

---

Published in Transactions on Machine Learning Research (10/2023)
0.0
0.5
7 0.9
Nomma mo 52D
mot PIGAM CUNGO
0.0
05
09
0.0
0.
50
0.5
0.9
(a) SHD
1.04
0.5-
FIR
0.2-
da
0.3
(b) FDR
0.5
0.6-
0.4-
05
05
05
09
(c) TPR
0.0
05
0.9
2.1-
25-
NNZ/P
20-
1--
3.0
63-
c0-
0.5
05
(d) NNZ/P
0.9
0012
0.0
-047
06
da
do
05
(e) FU for MN-CASTLE
Figure 16: The figure depicts results returned by additional monitored metrics for the considered (7, p)) settings. Here we use box plots, where we overlay the values of the metric, plotted as points.
32

---

Published in Transactions on Machine Learning Research (10/2023)
I=
0.9
1.0
0.5
Kendall-T
0.0
-0.5
-1.0
0.5
0.9
0.5
0.9
0.0
0.5
(a) Kendall-7
0.0
0.5
0.9
1.0
0.5
Spcarman ps
0.0
-0.5
0.0
0.5
0.9
0.0
0.5
0.9
0.0
0.5
0.9
(b) Spearman's rank correlation
0.0
0.5
0.9
1.0
0.8
nDCG
0.6
Model
0.2
MN-CASTLE
Random
0.0
0.0
0.5
0.9
0.0
0.5
0.9
0.0
0.5
0.9
(c) nDCG
Figure 17: The figure depicts box plots along with quartiles reference lines (dashed lines) for (a) for Kendall-r metric, Spearman's rank correlation, and (c) normalized discounted cumulative gain (nDCG) at 5. MN- CASTLE is given in blue while random baseline model in orange. For every dataset generated according to given (T) setting (i) sample 1 10 causal orderings 2 PL(O), where 0 is the estimated vector of scores; (ii) we draw x 10 random causal orderings ス ~ PL(0), where 8, ~ UO,, = .....N. Afterwards, we evaluate the three monitored metries by using the sampled causal orderings and - for both models. Thus, each box plot is made by 2 x 10a points.
33

---

Published in Transactions on Machine Learning Research (10/2023)
2.5
MI-CASTLE
Normalzed SHD
2.0
MSCASTIE
0.5
0.0
0.25
(a) SHD
11.8
01.6
FDR
0.6
0.2
0.25
0.75
(b) FDR
1.0
0.8
TPR
11.6
0.4
0.2
0.25
0.75
(c) TPR
CASTLE MSCASTIE E GOLEM
2.5
WNZIP
1.5
1.0-
0.25
0.75
(d) NNZ/P
0.00
0.00
0.00
0.25
0.75
(e) FU for MN-CASTLE
Figure 18: The figure depicts results returned by additional monitored metrics for different MN-DAG den- sities 8. Here we use box plots, where we overlay the values of the metrics, plotted as points.
34

---

Published in Transactions on Machine Learning Research (10/2023)
Noemalized SHD
10
15
20
N
(a) SHD
0.6
FDR
0.4
0.2
0.0
(b) FDR
7RR
01.6
0.4
0.2
10
15
20
N
(c) TPR
N
(d) NNZ/P
0.00
0.00
0.00
(e) FU for MN-CASTLE
Figure 19: The figure depicts results returned by additional monitored metrics for different number of nodes N. Here we use box plots, where we overlay the values of the metrics, plotted as points.
35

---

Published in Transactions on Machine Learning Research (10/2023)
1.0
Normllaxd SHD
0.65
0.2
01.0
No violation
Laplace
Undorm
Subsampling
(a) SHD
1.0
0.8
FDR
11.6
Laplace
Subsampling
(b) FDR
1.01
0.8
7PR
11.0-
0.6
0.2
Laplace
Uniform
Subsampling
(c) TPR
NNZ/P
1.0
0.8
violation
Laplace
Uniform
Subssmpling
(d) NNZ/P
0.5
0.6-
0.6
0.2
0.0
No violation
Laplace
Uniform
Subsampling
(e) FU
Figure 20: The figure depicts results returned by additional monitored metrics for different kind of generative model misspecification. Here we box plots, where we overlay the values of the metrics, plotted as points. To case the comparison, also provide the case without any violation.
36

---

Published in Transactions on Machine Learning Research (10/2023)
10-
0.5-
Kandall-r
0.-
-1.0-
0.25
0.5
0.75
15
Laplace
Uniform Violation
Subeampling
(a) Kendall-r
0.5-
Spaarman's &
0.0-
-0.5-
0.25
0.5
0.75
15
Laplace
Unitoem Violation
Suhsampling
5
(b) Spearman's rank correlation ps
1.00
0.8
+DCGaN
0.06
0.2
MN-CASTLE
Random
0.0
0.25
0.5
0.75
10
15
20
Laplace
Uniform
Subsompling
Violation
(c) nDCG@5
Figure 21: The figure depicts box plots for ja Kendall-r, (b) Spearman's rank correlation, and cc the normalized discounted cumulative gain (nDCG) at N. MN-CASTLE given in blue while a random basclino model in orange. Box plots on the left refer to the second experimental configuration, i.e., when we vary 8 while keeping fixed the values of the others parameters, as described in [Section 7.2. Box plots on the center concern the third experimental setting, where study the sensitivity of the cstimation accuracy w.r.t. the network size N. Box plots on the right relate the fourth experimental setting, where we study MN-CASTLE under generative model misspecification. Notice that in the left and right plots N = 5.
37

---

Published in Transactions on Machine Learning Research (10/2023)
Finally, Figure 21 depicts the results for Kendall-T statistics, the Spearman's rank correlation, and the nDCG at N related to causal ordering estimation. Here, we see that MN-CASTLE, performance grows along with 6, does not show dependence (in mean terms) on the size of the underlying MN-DAG, and does not decrease under generative model violation.
Appendix K
Data sources, Pre- Pre-processing and Inference Details
We obtained data on a weekly basis from January 1, 2018 to December 31, 2022. In detail, we downloaded Henry Hub natural gas futures prices (NG), WTI futures prices for crude oil (CO), New York Harbor No. Heating Oil futures prices (HO), and US natural gas storage (ST) data from the website of the US Energy Information Administration (EIA). The crack spread was calculated as the difference between HO and CO, with HO being converted to dollars per barrel. The deviation of storage from the norm (SDD) determined by comparing the ST value for given week to the average value for the same week over the previous five years. We also downloaded rig counts (RC) data from Baker Hughes, and extracted deviations from seasonal average values of gas consumption for cooling (CDD) and heating (HDD) environments from the National Occanic and Atmospheric Administration website. Finally, the economic uncertainty index (UT) was downloaded from the Federal Reserve Economic Data repository.
We have checked for non-stationarity in the time series using the ADF test at a 99% level of significance. We have found evidence of a unit root in all time series except CDD and HDD. Therefore, we have taken the first differences of the non-stationary time series. Finally, to account for differences in scale among the time series, we have standardized each time series.
After the preprocessing, we have obtained an estimate of S, by using the R package mvLSW ( Tayylor et al. 2019)). The wavelet transform has been performed using Daubechies wavelet with filter length equal to In addition, the smoothing of the periodogram has been performed using the rectangular kernel, also known Daniell window, of width 16. Finally the smoothed periodogram has been corrected for the bias by using the inverted autocorrelation wavelet ([Eckley & Nason 2005)) and regularized to ensure positive definiteness (/Schnabel & Eskow 1999)). At this point, we have obtained an estimate of Oj by means of inversion for the first 5 finer scales. Indeed, the values for the estimate S, were not negligible over these 5 scales.
In order to estimate the mode of the PI. distribution 40, we have proceeded follows. We have run the first step of the inference procedure have computed the estimate i by calculating the median over the previous estimates. At this point, we have obtained ="= argsort(-0). Indeed, since our model is probabilistie, we may have obtained different estimates of 0 at each run.
1000 times, and we have obtained n estimates O,, i E [n]. Then, we
Finally, to infer the causal coefficients we have used combination of kernels, specifically K K Periodic + K Linear X K 3/2 and we have set p 0.05 for hard-thresholding the partial coherence (sec Section 6).
38

---

Published in Transactions on Machine Learning Research (10/2023)
Appendix L
Complete MN-DAG of the Real-world Application
Scale
Scale
HOD
HPD
R8C
RC
ED0
CDD
NG
Scale
Scale
H00
HDD
RC
RC
CA0
0A0
NG
NG
C5
8
Scale
Storage deviation: SD
HDD
Uncortainty
Heating days norm deviation: HDD
Rig count: RC
Cooling days
CO0
Natural Gax: NG
NG
Crack sproad: C5
Figure 22: The figure depicts the multiscale DAG retrieved by MN-CASTLE, where nodes are sorted ac- cording to the estimate 0
39


Reviews:
[{"review_id": "3dwl33PHzt", "reviewer": "Reviewer_KzQ3", "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": "Strength: \n\n- The paper discusses both the generation and learning of the MN-DAG.\n- the problem is not carefully studied before, and hence a good novelty\n- proposed method is non-trivial, although most are based on components of some existing DAG approaches. \n\nWeakness:\n\n- No guarantees on the soundness or consistency of the learned graphs. \n\n\nSome comments/questions: \n1. time scale 2^j: does it mean there are two graphs for j = 1, 4 for j=2 etc?\n2. \"Causal inference methods\": typically these are called causal discovery methods to infer graphs, inference are typically limited to parameter and ATE as such.\n3. Section 3 first paragraph: seems like an exact repeat of a paragraph in the introduction.\n4. Proposition 4.3: How does one understand the spectral properties of X_T? is it the graph property?\n5. Figure 4: caption \n6. Does the policy gradient based learning method require a lot of data samples?\n7. It was not fully clear to me how the hyperparameters $\\mu$ and $\\tau$ are used in the experiments. Appendix H should be the main text of the paper on estimation, and does not discuss $\\mu$.  Are $\\mu$ given or estimated?", "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 proposes a graphical model, MN-DAG, to multiscale and nonstationary time series datasets. Generation of MN-DAG uses hyperparameters to model different time scale and nonstationarity among data. An variational Bayes method to infer the graph structure is provided, based on causal order and MN-DAG. Experimental results show promising performances of the proposed methods over a set of baselines.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Besides the question above, the organization of the paper should be more carefully thought. Before talking about Sampling an MN-DAG, one section should discuss the model itself (instead of in the introduction). The order of data generation and then learning also seems confusing.", "broader_impact_concerns": "n/a", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "5ksoebDC8m", "reviewer": "Reviewer_URns", "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- The model is ambitious and general in permitting both nonstationarity and the multiscale property\n- Extensive theoretical background is included\n- Real world experimental results demonstrate the effectiveness of the method \n\nWeakness\n- While the method contains some novel improvements (which improves from the previous submission), the generative model and learning procedure are based on off the shelf methods and frameworks so novel content is somewhat limited", "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 proposes MN-DAG, which extends structural causal models to include nonstationarity and a multiscale property (causal relations changing over time).\n\nA generative model is describe which involves (i) sampling the number of time scales from a binomial distribution, (ii) sampling a causal ordering according to the Placket-Luce distribution for each time scale, (iii) sampling and masking to generate weight matrices. Data may then be sampled from this model using the locally stationary wavelet process.\n\nThe authors then describe MN-CASTLE, an algorithm for learning MN-DAGs from such data using stochastic variational inference.\n\nExperiments show data sampled from the MN-DAGs have the expected properties, MN-CASTLE is more effective at learning from synthetic data simulated from using the described procedure when the multiscale property is high and MN-CASTLE is able to estimate multiscale causal structures in the US natural gas market that vary at the beginnings of the COVID-19 pandemic and Ukraine war.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "This submission represents an improvement both in terms of novelty (Lemma 4.2 and Proposition 4.3) and experimental results (overall improvement and the inclusion of the US natural gas market example which exhibits the multiscale property in real world data), which were my primary two concerns when reviewing an earlier version of this paper. \n\nI would still, however, recommend that in the introduction the authors replace 'structural equation models' with 'additive noise-based' or something similar to distinguish ANM / LiNGAM approaches from constraint-based and score-based structure learning methods, since the later two refer to structure learning and 'structural equation models' is typically used to refer to the framework more generally (not a specific learning procedure).", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "7m3UiEIqf1", "reviewer": "Reviewer_Luq6", "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": "Strength:\n- The paper is well written and clear.\n- The topic considered in this work is interesting and important because it is common to encounter non-stationary in real world.\n- The proposed model and estimation method seem to be sound to my knowledge.\n-  Experiments are conducted across different settings and the results are convincing.\n\nWeakness: My major concerns have been addressed by the resubmitted manuscript. There was some concern related to identifiability of the proposed model, which the paper discussed as part of future research directions in Section 8. ", "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 studies the problem of causal structure learning from time-series data. This paper introduces multiscale non-stationary directed acyclic graphs (MN-DAGs) which generalizes linear DAGs to the time-frequency domain. The authors  develop a probabilistic generative model that allows sampling MN-DAGs based on user-specified priors on the time-dependence and multiscale properties of the causal graph. A Bayesian method (called MN-CASTLE) for learning the MN-DAGs is further proposed by using stochastic variational inference. Experiment results on synthetic and real world data are provided.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "As mentioned, my major concerns have been addressed by the resubmitted manuscript.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: I apologize for the delay in this decision. The reviewers and I are in agreement that this paper should be accepted to TMLR, and I recommend accepting as-is.

The paper has improved dramatically from its original submission, and the reviewers remaining concerns about this version have been addressed by the authors. In particular, the inference methodology has been improved and the experimental section expanded and made more realistic. The updated paper is well written, of interest to the community, and describes novel and relevant work.
Decision: Accept (poster)