HiBaNG: Hierarchical Bayesian Nonparametric Granger Causal Discovery in Low-Data Regimes
Authors: He Zhao
Vassili Kitsios
Terence O'kane
Edwin V. Bonilla
Venue: Accepted by TMLR ()

Abstract: We present a principled probabilistic framework for discovering Granger causal relationships from multivariate time-series data in low-data regimes, where short sequences limit the applicability of modern deep learning approaches. While deep neural vector autoregressive (VAR) models perform well in high-data settings, they often struggle to generalize with limited samples and provide little insight into model uncertainty. To address these challenges, we introduce HiBaNG, a hierarchical Bayesian nonparametric framework for Granger causal discovery. HiBaNG places a hierarchical factorized prior over binary Granger causal graphs that encodes structured sparsity and enables interpretable, uncertainty-aware inference. We develop a tractable Gibbs sampling algorithm that exploits conjugacy and augmentation for scalable posterior estimation. Extensive experiments on synthetic, semi-synthetic, and real-world climate datasets demonstrate that HiBaNG consistently outperforms both classical and deep VAR baselines, achieving improved accuracy and calibrated uncertainty.

Published in Transactions on Machine Learning Research (02/2026)
HiBaNG: Hierarchical Bayesian Nonparametric Granger Causal Discovery in Low-Data Regimes
He Zhao
he. www.tm1..cii..ooooooo
CSTRO's Data61, Australia
& Department of Data Science and A1, Monash Universily, Australia
Vassili Kitsios
mussili.kitsio acsiioo.
CSIRO's Environment, Australia
& Laboratory for Turbulence Research in Aerospace and Combustion, Department of Mechanical and Aerospace Engineering, Monash University, Australia
Terence J. O'Kane
terence. okanee @csiro.ae
OSIRO's Ensironment, Australia & University of Tasmamia, Australia
Edwin V. Bonilla
adwin.bonille @csiro.co on
CSIRO's Dala61, Austrulia
Reviewed on OpenReview: https: /// openreview. net/ form?? id= e4 R03Y1RBB
//
Abstract
We present principled probabilistic framework for discovering Granger causal relationships from multivariate time-series data in low-data regimes, where short sequences limit the applicability of modern deep learning approachos. While deep neural vector autoregressive (VAR) models perform well in high-data settings, they often struggle to generalize with limited samples and provide little insight into model uncertainty. To address these challenges, we introduce HiBaNG, a hierarchical Bayesian nonparametric framework for Granger causal discovery. HiBaNG places a hierarchical factorized prior over binary Granger causal graphs that encodes structured sparsity and enables interpretable, uncertainty-aware inference. We develop a tractable Gibbs sampling algorithm that exploits conjugacy and augmentation for scalable posterior estimation. Extensive experiments on synthetic, semi-synthetic, and real-world climate datasets demonstrate that HiBaNG consistently outperforms both classical and deep VAR baselines, achieving improved accuracy and calibrated uncertainty.
1 Introduction
Multivariate time-scrics (MTS) data consist of obscrvations of multiple variables recorded at multiple timestamps and fundamental to a wide range of applications in economics, healthcare, climatology, and neuroscience. In these domains, uncovering causal relationships among time-series is often essential for understanding system dynamics and supporting decision-making,
In this paper, we focus on discovering such relationships from observational MTS data using Granger Causality (GC) Granger 1969 Lütkepohl] 2005 Shojaie & Fox 2022), which posits that one variable is causal for another past values provide statistically significant information for predicting the future of the latter, beyond what is contained in its own past. While GC does not necessarily imply a structural interventional causal relationship in the sense of, e.g., Pearl's do-caleulus (Pearl 2009) Rubin's potential outcomes (Tmbens) & Rubin 2015) it remains widely used for time-series causality The most common implementation of GC
Throughout this paper, causality refers specifically Granger causality.

---

Published in Transactions on Machine Learning Research (02/2026)
is via Vector Autoregressive (VAR) models [Lütkepohl 2005 which assume that each variable is function of the lagged values of other variables.
Bayesian VARs [Wozniak] 2016
Miranda-Agrippen & Ricco 2019) extend classical VARs by incorporating
prior beliefs, enabling uncertainty quantification and improved estimation under limited data. More recently,
deep parning-based VAR models (Montalto et al. 2015] Tank et al. 2018 Wang et al., 2018 Nauta et al. 2019; Khanna & Tan 2020 Wu et al. [000]] MMrriinkvvis & Vogt 2021 Gong et al., 2022 Fan et al.) 2023 have gained popularity for their ability to capture complex, nonlinear dynamics, provided ample data are available. In practice, however, many real-world applications do not have the luxury of abundant data, especially during the early stages of data collection. In these low-data regimes, where the number of samples is small relative to the number of variables, deep VARs often underperform. This is mainly due to the following: 1) When model complexity exceeds the informational content of the data, the learned function becomes unstable, noise-sensitive, and yiclds unreliable predictions with distorted uncertainty estimates. (Geffner et al. 2022 Annadani al., 2023 Deleu et al., 2023): 2) In practice, model selection and hyperparameter tuning are difficult. in causal discovery because the task is unsupervised: the target object (the true causal graph) is typically unobserved, and therefore cannot be used directly a validation target Biza et al. (2020). Morcover, although ground-truth graphs may cxist in synthetic benchmarks, using them for model tuning can lead to potentially biased comparisons, since it amounts to selecting models with privileged access to evaluation information that is unavailable in real applications Machlanski et al. 2023
To address these challenges, we propose a new hierarchical Bayesian VAR framework for Granger causal discovery, specifically designed for low-data regimes. Our method addresses the aforementioned issues by:
1. Introducing a hicrarchical factorized prior over binary Granger causal graphs, which is designed to encode structured sparsity, adaptive regularization, and low-complexity latent organization, helping stabilize causal discovery in low-data regimes where likelihood-basse evidence alone insufficient;
Decomposing Granger causality into discrete (graph structure) and continuous (causal strength) components, allowing the binary graph to constrain parameter estimation and improve generalization;
3. Leveraging Bayesian nonparametric techniques to integrate out latent factors and reduce the number of tunable hyperparameters.
Overall contribution: We develop HiBaNG, a hierarchical Bayesian nonparametric framework for Granger causal discovery in data-scarce settings, which integrates interpretable priors, principled uncertainty quan- tification, and enables tractable posterior inference via Gibbs sampling. Our extensive experiments on synthetic, semi-synthetic, and real-world climate datasets show that HiBaNG attains improved or competitive performance relative to both classical and deep VAR baselines.
2 Preliminaries
Consider collection of MTS data consisting of N time series or variables and T timestamps, stored the matrix X E NxT (Ø) (.1...,, RT), where Φ, E N consists of the samples/vall /values of the N variables at timestamp t E {1, fT). T). assumes that D, can be predicted from the mmax time lags (#;-11.... -max by learning coefficient matrix A E NNN for each lag E fl...., m }:
A VAR model for Granger causality (Lütkepohl 2005 Hyvärinen et al. 2010)
2 B + Et:
(1)
where et is an independent noise term. Conventionally, variable j (the parent) does not Granger-cause
variable i (the child) (i.j E (1....,N) N)) if and only if for all T, A7, = 0. For deterministic VARs, learning
can be done by minimizing a regression error: min ||2; Eimr A + X reg([AT) Tmmx)
(mmx)) ), where
reg( A Jr TTmax) ) is a sparsily-inducinng penalty (e.g., a group lasso penalty (Yuan & Lin 2006
Lozano et al.
2009 ): >N Ex II I12. Other alternative penalties can be found in Nicholson et al. 2217
We superscripts index matrices, rather than denote powers.
2

---

Published in Transactions on Machine Learning Research (02/2026)
another important line of research, especially in econometries and finance. A standard method may model with a multivariate normal distribution: Tmax 2; MM e Lrll Ar 2), where various priors can be imposed on [AT)E (e.g., -indutiigni Wiinsrr priors). Learning BVARs involves inferring the posterior of {AT] Tamx and Σ. In standard BVARs, one may need to "convert." " {AT) J
Bayesian VARs (BVARs) Litterman 1986)
priors) and E (c.g.,
Tman into GC graphs.
Deep VARs have recently become popular, as they use deep neural networks to model nonlinear dynamics across timestamps, which essentially generalize Eq. (1) using: zit fi (E Fmar A A' Ct where fi is typically implemented nonlinear neural network.
3
Method
In this section, we first give a complete picture of our proposed model and then discuss the motivations for our design choices at the end of the section. Specifically, we first present our BVAR model in a general form, which separates the coefficients into binary GC graphs and weight matrices; then we propose a new link function that helps build a hierarchical model on binary GC graphs in Section 3.1 subsequently, we provide details of the full model in Section 3.2 and its Bayesian inference algorithm in Section 3.3 Finally, in propertics and design choices in our model that allow us to tackle the challenges
Scction 3.4 discuss mentioned in the introduction.
Our overall model can be seen as a Bayesian VAR of the form:
It ~ MN
OA
(2)
where GT E 10, 1xN is the adjacency matrix for the binary GC graph of lag and @ denotes the Hadamard product. We further impose the following conjugate prior distributions (Miranda-Agrin Agrippino & Ricco 2019 A:: Vi, Gamma 1), AT, A (0, (+I,)=") and E: X Gamma (1, 1). ) = diag (A1,..., AN) where diag (A1....,AN returns matrix with its diagonal clements as AI,..., AN. In our model, the impact of variable j on i is modelled by two components: GT E {0,1) }} indicating whether there is a link bet and j and AY E indicating the weight of the link. If GT = 0, j does not impact i in lag 7 regardless of the value of A while if GT, = 1, A captures the influence from j to i. We also note that variable j does not Granger-cause cause variable (i.j E (I....,N)) and only if for all t, G O.
3.1 Generalized Bernoulli Poisson Link
Before describing the model in full detail, we propose a new link function named generalized Bernoulli Poisson link (GBPL) that thresholds a random Poisson variable m at V E (1,2.... } to obtain a binary variable b, one of the key building blocks in our approach.
Definition 3.1. (Generalized Bernoulli Poisson Link)
Poisson (1) b 1(m 2 V),
where 1(-) is function returning one if the condition is truc and zero otherwisc.
Property 3.2. Given 7 and V, one can marginalize out to get: b ~ Bernoulli (1 (1- 11 _ -77
V--
Remark. As b = if and only m < V, p(b 0) = EV- p(m = v), thus, p(b = = 1) 11 DV- p(m = v). Moreover, as E(b) = p(b 1), larger leads lower expected probability of b being one, under the same 7 Property 3.3. Given b, the conditional posterior of conditioned on b, 7, and V is in closed form:
TPoissonv (7),
b 0
m
Categoricaly ([... f(.,)).... otherwise
where TPoissony (7) is the Poisson distribution with parameter y left-truncated at V (i.e., samples from that
-AA="
Poisson distribution are greater than or equal to V) and [(u,7) = ΣV -4
Ager is the normalized Poisson
3

---

Published in Transactions on Machine Learning Research (02/2026)
probability mass function. To sample from truncated Poisson distribution, common approaches are rejection sampling (Geyer) inverse transform sampling with inverse CDF (NumPyro). Although either approach is efficient when V is large, as shown later, V in our case takes a small number, leading to relatively efficient sampling algorithm.
Remark. If b > O, m 2/V V almost surely (a.s.) and one can sample m from the truncated Poisson distribution at V efficiently by computing the inverse Poisson cumulative distribution function GGiles 2016). If b = 0, m E {O..... V- 1} is sampled from the categorical distribution with normalized Poisson probability masses. The closed-form conditional posterior contributes to the development of an efficient algorithm of our model.
Property 3.4. When V is set to 1, GBPL reduces to the link function proposed in Zhou 2015
3.2 Poisson Factorized Baange-Caussl Graph
Now we introduce our Bayesian construction on binary GC graphs with GBPL. To assist clarity, we discuss our method with only one lag, i.e. Tmax 1, temporarily omitting the notation of lag T, and introduce the extension to multiple lags later.
The general idea is that we assume binary GC graph G E 11NNN is sample of a probabilistic factorization model with K latent factors: G~ p(O") where O R Nxk each entry of which 0i indicates the weight of the Athh factor for variable i of being a child in GC relation and & E R each entry of which 05,s indicates the weight of the 4th factor for variable j of being a parent in a GC relation. In this way,
whether j Granger-causes i depends on their interactions with all the K factors: Gi, P (E (E 0c.005,.)
Conditioned on θ and 2 we have: p(G 2) = II I[ p(G,,, , 2), meaning that the links in G can be gencrated independently.
With the help of GBPL, we propose to impose the following hierarchical Bayesian prior p(G):
Tk Gamma (1/K,1/0) 0ix ~ Gamma (a;, 1/dx), 03.h Gamma (bj, 1/ex),
M.j Poisson
Σ TklLLOs
Gas= = 1 2 v),
(3)
where the first and second parameters of the gamma distribution are the shape and scale parameters, respectively; noninformative gamma priors Gamma (1, are used for a, by, de, cke and c.
In the above model, variable Fk is introduced to capture the global popularity of the 4thh factor (Yang & Leskoveo 2012 2014 Zhou 2015). Intuitively, we can understand the vectors 0g,, and dj., as the embeddings of variables and j, respectively. Whether i is child of j is determined by the inner product between their embeddings. Morcover, each dimension of the emboddings weights differently and r& indicates the weight of the hthh dimension.
Mathematically, the construction on θ 4, and r can be viewed as the truncated version of a gamma process (Ferguson 1973 Wolpert et al.) 2011 Zhou 2015) on a product space R, x Ω: & ~ IP(G,,, 1/c), where Q is a complete separable metric space, c is the concentration parameter, Ba,s is a finite and continuous base measure over Ω. The corresponding Lévy mcasure vdddddd)) re- dr (dDdd)) In our case, draw from Ba,b is pair (0.,,,.)) where 0,, = [0I,....., , ON and Φ = [ol.,,..., &N,A A draw from the gamma process is discrete distribution with countably infinite atoms from the base measure: 6 E r. and F is the weight of the 8th atom. Although there are infinite atoms, the number of atoms with ra greater than p E R4 follows Poisson[[ - and the expectation of Poisson decreases when p increases. In other words, the number of atoms that have relatively large weights will be finite and small, thus, a gamma process based model has an inherent shrinkage mechanism. In our case, if we set the maximun number of latent factors K (i.e., the truncation level) large enough, the model will automatically learn the number of active factors.
In our model, Bayesian nonparametrics provides a principled way to build models whose complexity adapts to the data, rather than being fixed ahead of time. Even if K is set to a large value, the model will not use all of the latent factors; the prior acts as complexity penalty, preferring simpler models unless there
4

---

Published in Transactions on Machine Learning Research (02/2026)
HIBaNG
Data
Gibbs Sampling
Multvrrtte Time Series
WWW M mmm
Hierarchical Factorized Prior
mm
Posterior over
Granger Causal Graphs
Figure 1: HiBaNG models the given MTS data via set of binary Granger causal graphs G" as latent variables with a hierarchical factorized prior.
is strong evidence to the contrary. In model, K is quite different from hyperparameters in parametric models such as the number of layers of neural network. K is truncation level that tells the model the maximum number of latent factors it can potentially use, while the actual number it uses is determined by the data. In addition, the model will self-regularize if it uses more latent factors than necessary.
Finally, we refer to the model in Equations (2) and (3) Poisson Factorized Granger- Causal Graph (PFGCG) and denote G ~ PV)CG(V). In the case of multiple lags, we summarize the model
GT (V) xt xx~ xM~
Σ O O G' O E
(4)
where we have a separate generative process for the GC graph at each lag T. The overall model is shown in Figure 1
3.3 Inference via Gibbs Sampling
Here we introduce how to estimate the posterior over the parameters of the above model using Gibbs sampling (Casella & George 1992), which adheres to the detailed balance condition (Gilks et al. 1995). a fundamental property of Markov Chain Monte Carlo (MCMC) methods that guarantees the Markov chain converges to the desired posterior distribution as its stationary distribution. To further enhance efficiency, our method employs hicrarchical prior structure with conjugacy properties, with the help of several augmentation
techniques between Poisson and gamma distributions (Zhou et al., 2012 Zhou 2015 These conjugate priors ensure that all conditional distributions are analytically well-defined and computationally tractable, simplifying the sampling process. The conjugate structure not only improves computational efficiency but also contributes to the stability Gibbs sampling, enabling the algorithm to effectively explore the posterior distribution even in challenging settings. Here we highlight the sampling G" and leave the other details in Scction
An entry G in G7 is involved in the generative process of data in Eq. (2) and has Bernoulli prior
according to Eq. (3 Therefore, by denoting p(GT,, = 0)-) = s the other variables) ), we can derive:
7.0 and p(G 1|-) = s (- stands for all
Th (gT,))
r.l - 3 227,, (1 -700)
8
and 8
(5)
where:
Σ
UJ= Σ
1t-7 W Σ Z4 Pzze-rs
t=l
N
Timax N Σ
List Σ
GT
(6)
-
5

---

Published in Transactions on Machine Learning Research (02/2026)
With the above conditional posterior, one can After each sample, we only need to update Wi, and the other statistics can be updated after all the entries are sampled.
We can then sample G ~ Bernoulli (
/ 1. +
+
sample the entries of G" one by one using Eq.
5
Computational complexity: With the above we have that, in one Gibbs sampling iteration, the complexity of sampling GO is O(N2). As can be seen in Algorithm [ in the appendix, the whole complexity of each Gibbs sampling iteration of our model is (N + K Tmax + T where N, T, V, K, Tmax are the number of variables, the number of timestamps, the truncation level, the number of factors, and the maximum mumber of lags, respectively.
3.4 Properties, Design Choices, and Practical Consequences
We refer our method as HiBaNG (for hierarchical Bayesian nonparametric Granger causal discovery), which uses the PFGCG model and the Gibbs sampling algorithm described in Sections 3.2 and 3.3 respectively. Here we discuss some of HiBaNG's properties, its underlying design choices and their practical consequences.
Separate A and GO structures modeling VAR coefficients: In our approach, A O G represents coefficients, which are intrinsically sparse as G is binary. Thus, no sparsity-induciig penalties Nichollon et al. 2017 Ahelegbey et al. 2016 Ghosh et al. 2018 Billio et al. 2019) or post-hoc heuristics (Nauta et al. 2019: Marcinkevičs Vogt 2021) are needed. Morcover, modeling binary G directly enables us to learn the probability of causal link between two variables directly and one can quantify uncertainty straightforwardlly
Independence of Granger graphs across lags: In our model, assume that the prior over Granger causal (GC) graphs factorizes across lags as shown in Eq. (2). Our primary motivation for this choice is to avoid imposing strong parametric or structural assumptions about temporal coherence of causal effects across lags, which may not be universally valid. While many dynamical systems exhibit smoothly decaying or persistent effects, others such as in complex observational domains such as climate or socio-economic systems can exhibit lag-specific interactions that do not conform to monotonic or smooth temporal patterns. An independent prior across lags therefore provides more agnostic and flexible inductive bias. Importantly, although the priors are factorized, the posterior is not. The lag-specific graphs are coupled through the shared likelihood: ( Ymur X) x p( (G(T) I p(G(T)), which couples all lagged parent sets through the autoregressive structure, inducing posterior dependencies across lags when supported by the data. Thus, temporal coherence is learned rather than imposed, which is especially important in low-data regimes where strong priors can dominate inference. Although cross-lag sharing can be introduced by tying parameters or latent variables, this adds modeling commitments and weakly identified hyperparameters that can increase sensitivity to prior choices. We therefore use a factorized prior across lags, allowing coherence emerge through the joint likelihood without additional tuning.
Linearity and potential extension to nonlinear models: Our model adopts a línear Gaussian VAR likelihood, which emphasizes interpretability, principled uncertainty quantification, and robustness in low-data regimes. In principle, the framework can be extended to a nonlincar likelihood of the form Φ; = 1 (1 Tmax [GT))rr) + where f is a parameterized nonlinear function (e.g,, neural network). However, such an extension would generally break the conditional conjugacy structure of the proposed Bayesian model. As a result, inference would require hybrid algorithms (e.g., combining MCMC with variational or gradient- based methods), substantially increasing computational complexity and reducing robustness. We note that a linear VAR may not be a faithful generative model for many real-world systems, and that fitting a lincar model to data generated by nonlinear dynamics can introduce bias. Nevertheless, one may consider using linear model even when the underlying system is nonlinear as pragmatic approximation. In such cases, linear Granger- Granger-based methods may still recover meaningful aspects of the causal structure when dependencies are approximately linear in expectation, or when nonlinear effects give rise to detectable linear predictive improvements. Moreover, in low-data regimes, linear models can offer more stable inference and belter-calibrated uncertainty than highly flexible nonlinear alternatives.
Suitability to low-data scttings: Our hicrarchical Bayesian model is particularly suited to low-data settings because it allows for information sharing across related parameters through structured prior. Specifically, our prior on G takes factorized form, allowing Oi or dj. to capture the specific information
6

---

Published in Transactions on Machine Learning Research (02/2026)
for individual variable in terms of factor k. At deeper level in the hicrarchy, 0a, is inffuenced by two higher-level components: dy (capturing variable-specific traits) and de (capturing factor-specific traits). This setup allows the model to partially pool information to learn robust estimates for 02 even when direct data is sparse by borrowing strength from related variables (via dk) and related factors (via as). As a result, the model remains expressive while enhancing generalization in low-data regimes.
Need for the GBPL link function: Instead of modeling binary GC graphs directly with Bernoulli distributions, we introduce GBPL as a link function that connects Bernoulli variables to underlying Poisson distributions (Zhou 2015). This transformation enables the use of hierarchical Poisson-Gamin Gamma constructions, akin to those in Poisson matrix factorization models |Canny) 2004; Zhou et al. 2012 Gopalan et al., 2014). By leveraging the Poisson-gamma conjugacy, we gain access to a broader and more flexible set of tools for hierarchical Bayesian modeling and inference, which are difficult to apply directly to Bernoulli likelihoods.
Given the specification of the prior distributions, one can see that E Ek=IFO0,, B((),)) = 1 and B(rx)) = 1/K. Therefore, according to Property 3.2 of GBPL, a priori, the expected sparsity of G is N² (1- -I- ) Note that V E . is a hyperparameter that incorporates our prior belief of the graph sparsity, i.c., larger value of V means that we encourage the model to learn sparser graphs, while the final sparsity will be determined by the model to fit the observational data. Specifically, we compute the value of 1 - -' e- of with V = 21, 3, 4} as 0.2642, 0.0883, 0.0190, 0.0037, respectively. When V V= 3, it means that less than 2% of the node pairs are expected to be connected in a graph. When V = 4, the sparsity is less than 0.5%, which can be overly sparse. Empiricalyy, we observe that when V ==4, the sampled GC graphs from the posterior nearly have zero links, thus, we set V e {1,2, 3} in practice. According to Property 3.4, the link function proposed in Zhou (2015) is special case of GBPL (when = 1). The model of Zhou 2015 cares less about the sparsity of a graph as it is given in the data. However, the expected sparsity of /Zhou (2015) is 0.2642 (V 1), which is too dense for many GC discovery problems. Although the expected sparsity takes finite values, the posterior graphs are not pinned to the prior grid, as data often overwhelms the prior.
as E(Or,s))
Fewer hyperparametors: In low-data scenarios without ground-truths, selecting a model from a large set of hyperparameters challenging. Our method has, by construction, a small number of hyperparameters integrates out the intermediate variables in its hierarchical Bayesian structure. Furthermore, it automatically learns K from data via Bayesian nonparametrics. Indeed, the main hyperparameter of our model is V that only takes a small number of discrete values, while the regularization weights in other VARs are usually continuous parameters in an infinite range.
it
4 Related Work
Here we focus on the following lines of related work in the machine learning literature and refer readers to surveys such as (hhoaaie& & Fox (2022) Assaad et al. (2022); ; Gong et al. (2023) for a more comprehensive overview.
Bayesian VARs: have been widely used in econometrics and statistics Breitung & Swanson 2002 George
et al.
2008
Fox et al.] 2011 Nakajina & West 2013 Ahelegbey et al. 2016 Ghosh et al. 2018 Billio
ct al.
2019
Ghosh et al. 2021). For comprehensive reviews, see Wožniak (2016); Miranda- Agrippino &
Ricco
2019 Unlike most. existing methods that focus on modeling real-valued VAR coefficients, our method
models the binary GC graphs directly via a novel construction. Since BVARs are studicd across multiple
disciplines, comprehensive comparisons to the latest deep VARs in the same settings and dalasets have not.
been carefully studied before.
Deep VARs: have recently become popular in the machine learning community Moontalto et al. 2015 Tank
el al. 2018 Wang el al. 2018 Nauta et al.) 2019 Khanna & Tan 2020 Wu et al. 2020 Marcinkevids & Vogi, 2021
Gong et al. 2022 Fan et al. 2023), where different neural network architectures or learning mechanisms
have been explored Tank et al. 2018 Nauta et al.] 2019 Khanna & Tan 2020 Marcinkovićs & Vogt 2021
Bussmann et al., 2021 Fan et al.) 2023 Zhou et al. 2024). More recently, Bayesian deep VARs have been
proposed, such as ACD (Löwe et. al. 2022 RHINO 0 (Gong et al.] 2022). Dyn-GFN (Tong et al. 2022
7

---

Published in Transactions on Machine Learning Research (02/2026)
and MCD (Varambally et al.,
2024 Although these methods can also be considered Bayesian approaches,
many of them have different methodologies and focuses from ours. For example, ACD (Lowe et al.) 2022) focuses on discovering causal relations across samples with different underlying causal graphs but shared dynamics. RHINO [Gong et al. 2022) extends VAR by modcling instantancous causal relations Runge 2020 Pamfil el al., 2020) and introducing bistory-dependent noise, which we do not consider in this paper.
Dyn-GFN (Tong et al. (022) is a Bayesian approach based on GFlowNets (Bengio et al.
2023) focusing on
discovering causal graphs varying with time. There are also recently proposed deep VAR variants focusing on causal discovery MTS data with missing values (Cheng et al. 2023 2024a). Wu et al.
(0024) studies a
different problem from ours, where the data are event sequences and each event has an occurrence timestamp and & type. While Deep VAR models focus on expressive nonlinear dynamics and high-capacity predictors, HiBaNG adopls a different objective: learning an interpretable and uncertainly-aware posterior over Granger causal graphs using a Bayesian framework that is well-suited to low-data regimes.
NoR methods: causal discovery on MTS data is not the focus of our paper. To capture instantancous
causal effects that are not modeled by VARs and GC, there functional causal models such as in Hyvärinen
et al.
(0100); Peters et al.
(2013); Pamfil et al. (2020) and methods based dynamic Bayesian networks
(BBNs) Dean & Kanazawa 1989: Murphy 2002) or structured VAR models in econometrics (Swanson &
Granger
1997
Demiralp & Hoover 2003). For DBNs, we refer readers to surveys such as Mihajlovic &
Petkovic
2001
Shiguihara et al. (2021). ). Moreover, there are also constraint-based approaches that extend
the PC hm (Spirtes et al. 2000 to model time-series data (Runge 2018 Runge et al. 2019 Runge
2020
Huang et al.
2020).
Causal discovery for non-time-series data: (e.g., i.i.d.) (Glymour et al. 2019) is another area. with
different focus. Here we consider Bayesian methods Lorch et al.) 2022 Cundy et al. 2021 Geffner et al. 2022
Bonilla et al.) 2024 Thompson et al. 2024 2025) that model binary causal graphs as loosely related work to ours among the rich literature. These methods usually leverage gradient-based Bayesian inference
algorithms such as variational inference and use reparamoterization techniques (Maddison et al., 2016 Jang
et al.] 2016 to relax the optimization over binary causal graphs to a continuous one, while ours models binary graphs directly. Moreover, they require the discovered graphs to be directed acyelic graphs. We believe that extending their methods to time-series data with multiple lags is nontrivial.
Other related works: Poisson factor analysis is generic Bayesian framework used in various areas such as graph learning (Zhou, 2015 Zhao et al. 2017p, topic modeling (Zhao et al.) 2018b a), and dynamical data modeling (Schein et al. 2019 To the best of our knowledge, it has not been adapted for Granger causal discovery. The closest works to ours include Kalantari et al. (2018); Kalantari & Zhou (2021). which learn a set of latent. factors from time-series data as well as a binary graph between them. There are several key differences with our method: 1) ours is tailored to Granger causal discovery, which learns binary graphs of time-series variables instead of latent factors; 2) ours considers multiple lags, while they only consider one;
and 3) ours is based on the proposed GBPL, while they use the original Bernoulli Poisson Link (Zhou 2015 which is less applicable to our problem.
5 Experiments
Here we evaluate our method using a set of synthetic and semi-synthetic datasets well as a real application involving the analysis of climate data. As mentioned in Section 3.4 we refer to our method as HiBaNG. In addition to the content of the main paper, we show comprehensive analyses of the empirical computational performance and convergence of HiBaNG in Sections D and of the appendix, respectively.
5.1 Experimental Settings
We use 10,000 as the maximum number of Gibbs sampling iterations, where the first 5,000 are burn-in iterations. We then collect samples from the conditional posteriors of the graphs every 10 iterations which are stored in Y E RR.N X x (H = 500 the number of collections). The Bernoulli conditional postcrior probability of a link between i and j at lag 7 in collection h E (1,...,FH H) } is computed by Eq. (5)
HiBaNG:
3 As shown in Section
D our method converges in fewer iterations.
8

---

Published in Transactions on Machine Learning Research (02/2026)
as /i + Given the collections, to compare with other methods, we compute the averaged probability of the discovered GC graph by mcan(max( dim 'T'), dim = 'h') (Marcinkevičs Vogt 2021). As HiBaNG has an intrinsic shrinkage mechanism on K, we set K 50, which is empirically large enough for our experiments. The only hyperparameter that we need to tune is V, which we vary in (1,
Baselines: As ours is VAR approach for GC, we mainly include baselines that are also based on the VAR framework. 1) We compare with the widely used VAR with F-tests for Granger causality and the Benjamini- Hochberg procedure (Benjamini & Hochberg 1995) for controlling the false discovery rate (FDR) (at q = 0.05), denoted as VAR (FBH) and implemented in the statsmodels library (Seabold & Statsmodels 2010). 2) For Bayesian methods, we compare with two classic approaches but with different prior distributions: BVAR with priors on the coefficients, named BVAR(d), i.c., (fA") J.mex, ; Σ) x |Σ] N+1
(Litterman 1986 Miranda- Agrippiio & Ricco 2019), whose posterior has an analytical form. BVAR with conjugate
buffuse/noninformaen
priors on the coeficients, named BVAR(c): vi,- ~ Gamma (L, 1), Ar, ~ N (0. (+) and for E: X This is equivalent to an ablation of our model without (GT),mm
Gamma (1, Σ diag for which we use Gibbs sampling with the same settings as ours. For deep/ /nep/neural VARs, we compare with a method with component-wise statistical recurrent units (SRU) (Oliva et al. 2017) and its improved version (economy SRU, eSRU) (Khanna & Tan 2020) with sample-efiiiint architectures. The important hyperparameters of SRU and eSRU are the strengths (11, 122, (aa) of three regularization terms. We also compare with GVAR (Marcinkevičs & Vogt 2021), which uses self-explaining neural networks (Alvarez Melis & Jaakkola, 2018 and converts the weights in the neural networks into binary GC graphs with a heuristic, stability-based procedure. As a state-of-the the-rtt method, Jacobian Regularizer- based Neural Granger Causality (JRNGC) (Zhou et al. 2024) is also included in the comparison. For noR methods, the most recent one PCMCI+ [Runge 2020) is also compared. For the baselines, we either use their original settings or follow those in Marcinkevićs & Vogt (2021), shown in Table 3 For all the compared methods, we set Tiax = 5
unless otherwise specified.
Evaluation metrics:
Following Khanna & Tan 2020 Marcinkevičs & Vogt (2021); ; which aggregate
graphs multiple lags into one, we use four metries to compare the discovered GC graph of method on a dataset with the ground-truthe -trut graph. For all the baseline methods, we compute the score of discovered GC graph from their learned VAR coefficients. For our method, the score of a GC graph is the mean of the Bernoulli posterior. We report the areas under receiver operating characteristic (AUROC) and precision-recall (AUPRC) curves by comparing the score of a discovered GC graph to the ground-truth graph. Moreover, as mentioned before, VAR (FBH) and GVAR use specific post-hoc processes to convert coefficients to binary GC graphs; thus, we also report the structural Hamming distance (SHD) between the discovered binary GC graph and the ground-truth one. Note that unlike AUCROC and AUPRC, SHD is biased toward the sparsity of the ground-truth graph; for example, for a sparse ground-truth graph, method that always predicts no links achieves a low SHD. To measure the predictive uncertainty of the Granger-causal graphs discovered by different approaches, we report the calibration error (CE) (Guo et al. 2017). which is a widely used metric for model uncertainty and confidence (Liu et al. 2020 Kumar ct al. 2019). CE examines the difference between the model's probability and the true probability given the model's 's output; its definition is shown in Definition 2.1 of Kumar et al. (2019). We treat the causal discovery task with N variables as a binary classification problem with N² samples (i.e., predicting a Granger-causs variables) and then compute CE accordingly using the method in Kumar et al. (2019). For AUCROC and AUPRC, higher values indicate better performance, and for SHD and CE, lower values are better. For all the numerical results, we run our method and the baselines with 5 different random seeds and report the mean and standard deviations.
link between pair of
Model selection and parameter tuning: Importantly, we note that our task is discovering GC graphs from data without training with ground-truth, which is an unsupervised problem (e.g., akin unsupervised clustering). As no ground-truth is given for training, we use forecasting performance, mean square error (MSE), on held-out temporal segments (last 20% of data (Gong et al 2022 as a proxy objective for selecting model parameters for VAR-based methods, including ours. This is consistent with how Granger causality is commonly defined based on predictive influence. We believe this constitutes a fair comparison across all VAR-based methods, as the same procedure is used consistently. For each method, we select the parameters
9

---

Published in Transactions on Machine Learning Research (02/2026)
of a method that give the best MSE. Our model selection is different from that of GVAR (Marcinkevičs & Vogt 2021), where the best model is selected by comparing with the ground-truth graphs and report the best achievable performance. PCMCI+ relies on conditional independence testing to infer a causal graph. Its key hyperparameter is the p-value threshold. Unlike VAR methods, PCMCI+ does not perform forecasting, and thus MSE cannot be used as a model selection criterion. Therefore, we used the slandard default setting of 0.05. The parameter space where we search for each method is shown in Table B of the appendix and the performance of our method under different parameter configurations is shown in Table | of the appendix.
Setting V: As stated in Section 3.4 the expected sparsity of G in the proposed prior distribution is N2 1 --- -V--- o) - where V controls the sparsity of the graphs. We compute the value of 1 - -YY- e- with V E [1,2, 4] as 0.2642, 0.0883, 0.0190, 0.0037, respectively. Thus, setting V=
of = 3 already induces significant sparsity, which matches the empirical needs of our experiments. We see that when V = 4, the sampled GC graphs from the posterior nearly have zero links, thus, we set V E (1, in practice.
Setting K: We believe that K = 50 is sufficiently large in our experimental settings. This can be seen on the RHS of Figure 2 which shows rx for all lags and all K factors as a K x mmax matrix (K = 50 and Tiax = 5). Recall that rhe indicates the weight of latent factor k at lag т. If we look at one column of the figure that shows the weights of the 50 factors for one lag, we can see that only a few entries have large values, meaning that only a few factors are active among K = 50. We had similar vations in other datasets as well.
5.2 Synthetic and Semi-synthetic Datasets
We report the results on toy synthetic datasets in Section B For more comprehensive quantitative comparisons, we conduct our experiments on three widely used benchmark datasots, detailed as follows.
Lorenz 96 [Lorenz 1996) is a standard benchmark synthetic MTS dataset for GC, which is generated from
the following nonlinear differential equations: dzis = (Zi+1.t 22) Ii,t + F, for 1 isis S N, where F is a constant that models the magnitude of the external forcing. The system dynamics become increasingly chaotic for higher values of F (Karimi & Paul 2010). We set N = 40, F == 40, and T= = (100, 500), which mimic noisy observations with reasonably large numbers of variables but few observations.
dt
Following Marcinkevičs & Vogt (2021), we evaluate the methods on another synthetic dataset generated by the Lotka- Volterra model (Bacaër & Bacaër 2011), where we use N = 40 and T = (200, 500). For the other parameters of the Lotka-Volterra model, we use the same settings as in Marcinkevičs & Vogt 2021
Finally, we consider the FMRI dataset with realistic simulations of blood-oxygen-lever dependent (BOLD) time series (Smith et al.] 2011). In the dataset, multivariate BOLD time series are generated from a known directed network structure via a biophysically grounded forward model that links neural activity to hemodynamic responses. Specifically, the simulator uses Dynamic Causal Modeling (DCM) style architectures (Friston) et al. 2003) to produce realistic BOLD signals from a predefined connectivity matrix. The corresponding directed adjacency matrix used in the gencration process is treated the ground truth when evaluating causal discovery methods. Following Khanna & Tan (2020 Marcinkevičs & Vogt (0021). we use 5 replicates from the simulation no. 3 of the original dataset, where N 15 and 7 = 200 are pre-specified as standard settings. We introduce straightforward metric to measure the severity of low-data conditions in a dataset: B = N/T, meaning that larger & indicates a dataset has a larger number of variables with a smaller number of observations. We show the value of / for each dataset in Table 1 We notice that Cheng et al. (2024b) recently introduced a few new benchmark datasets with many timestamps (c.g., from 8,000 to 50,000). Our focus is on low-data regimes with fewer timestamps, and their ground-truth graphs are undirected (i.e., the adjacency matrices are symmetric), while our method discovers directed graphs; therefore, these datasets are less applicable to our problem.
We present the results for AUCROC and AUPRC in Table I Overall, our proposed method, HiBaNG, demonstrates superior performance across most datasets and metrics. A notable observation is the enhanced performance of HiBaNG as the parameter & increases. This trend underscores the robustness of our model in severe low-data situations, where traditional methods often struggle. For the FMRI dataset, where the severity of low-data conditions is reduced, HiBaNG ranks second in both AUCROC and AUPRC. This result aligns
10

---

Published in Transactions on Machine Learning Research (02/2026)
Table 1: AUCROC and AUPRC results on semi-syuthetic datasets. VAR (FBH) and BVAR(d) failed to learn when T 100 on Lorenz 96 and T 200 on Notka-Volterra, Best and second best results are highlighted in boldface and underline texts, respectively. / N/T indicates the ratio of the number of variables to the number of observations.
Lorenz 96
Loreia 96
Lotka- Volterra
Lotka-Vollers Volteern
FMRI
T= 100, 3 0.4
7 500, B = 0.08
T= 200, / == 0:2
T= 500, 8 0.08
8 0.076
AUCROCT
AUPRCI
AUCROCT
AUPRC*
AUCROCT
AUPRCI
AUCROCT
AUPRCT
AUCROCT
AUPRC1
RVARS:)
0.47+0.02
0.09+0001
0.73+0002
0.43+0.02
0.50+0001
0.50+0.011
0.78-0.04
0.52+0.05
0.66+0.08
0442-0.12
BVAR(()
0.73+0003
0:43+0.02
0.67-0.01
0.23-+0.03
0.68+0.06
0:40-00066
VAR (FBR)
0.72+0.01
0.30+00.3
0.68-0.03
0.18+0.01
0.60+0.04
0.32-0.02
POMCI
0.2+00022
0:12-0:01
08240002
0.30+0.02
0.72+003
0.44+0.02
0.78-0.01
0.47+0.01
0.80+0004
0.67+0.07
SRU
0.5310.01
0.12:10.11
0.522:.03
0.57100.4
0.55::0.02
0.5010.01
0.81=0002
0.32.0002
0.6610.02
0.32
<SRU
0.3420.00
0.12
0.84:00011
0.6340.00
0.6410.03
0.5810.01
0.87=0.03
0.3611001
0.7240.01
0.47=0.01
GVAR
0.5740.01
0.15
0.5320001
[.6340.00
0.68:0.02
0.62400001
0.81=0.53
0.61440.04
0.7240.02
0.57=0006
JBNGC
0.5840.07
0.14
0.79:0.01
0.55200.3
0.6840.02
0.3240000
0.87=0002
0.61
0.6940.03
0.45=0004
HiBeNG
0.71:0008
0.35=00M
0.86:20.02
0.68:0.02
0.73:0001
0.78:0.01
0.81=0.88
0.66:0003
0.7320.03
0.60:0000
with expectations, as richer datasets reduce the relative advantage of our Bayesian approach. Comparing HiBaNG with BVAR(c), which essentially represents HiBaNG without the integration of binary GC graphs reveals that HiBaNG consistently outperforms BVAR(c). This comparison highlights the significance of incorporating binary GC graphs into the Bayesian framework, facilitating clearer and more interpretable understanding of causal links while enhancing prodictive accuracy.
The results for SHD and CE are displayed in Table 2 For SHD, we include methods capable of converting their coefficients into sparse graphs, while for CE, we concentrate on the top-performing methods based on AUCROC and AUPRC scores. It is important to note that SHD is inffuenced by the method's approach to generating binary graphs, which may rotuee bias depending on the sparsity of the ground-truth graphs. Our method's capability to directly sample binary graphs without relying on arbitrary thresholds provides distinct advantage. Regarding CE, our method achicves the lowest crror rates across almost all datasets. This indicates that HiBaNG's predictive confidence is highly aligned with its accuracy, result of its inherent uncertainty-aware design. Such aligmment is crucial in roal-world applications, where understanding the reliability of causal predictions can inform better decision-making.
Table 2: SHD and CE on semi-synthetic datasets.
Lecena
Lotka Volterr
FMRI
SHDL
VAR (FBI)
58 2.4
7420011044
28.81.
PCMET+
406.0:15.0
411.0:5.0
491.1120.0
423.0120.0
70.00:2.00
Figure 2: Qualitative analysis of HiBaNG on Lotka-Volterra, Left: Ground-truth GC graph, middle: Bernoulli posterior mean of the discov- K, show- ered GC graphs, right: Matrix of {r [rikk=
GVAR
$89.61220.6
127. 4L77.8
279001106.1
828124.1
71662218
IBDaNG
117.7-333
71.814.0
67.012.4
45.01-4.4
24.2.0.7
CEL
BVAR (d)
0.3020.01
0.112:001
0.1120002
ing that K = 50 is sufficiently large in our experi- ments (see Section 5.1 for more details). Each rect- angle indicates a value of a matrix and brighter
POMEI
0:200000111
0:27
0.2700011
0:28810111
03120002
GVAR
007=00011
0.36+00.11
0.08000011
0.3080011
0.190001
JRNGC 0.01110.1
0.0010.011
0.07100011
0.13110.01
0.09100022
ILBANG
0.11210011
0.08.10.01
0.05.10.01
0.0420001
0.07.0.01
colors indicates larger values.
we compare the Bernoulli posterior mcan of HiBaNG with the ground-truth graph on Lotka- Volterra. We can see that the posterior mean discovered by our method is well aligued with the ground-truth graph, where brighter rectangles indicate a higher probability of a GC link between two variables. Finally, recall that rf Eq. (3) models the weight of latent factor at lag т. We plot [rī] 1=1,,, matrix. The matrix is quite sparse: only a few entries have large values, and only a few factors are active among K 50. This demonstrates the shrinkage mechanism on K.
In Figure 2
K. as a K x Tmax
11

---

Published in Transactions on Machine Learning Research (02/2020)
M0-1-,
NAO (+
@NAD (1)
PNA PSSA (1-)
BURM (-1)
MME (1-2
MMME (1-1)
SCANNS 1h-1)
MEET FRMM2 (33)
MMB (5-5)
PSA1 (1))
8ID0
AR (i-)
RMMM (A))
BAM (+-1)
NNAO- MEE
DA
NNAD (1
PNA ()
PSA1 ()
RMMM ()
AA
1000
RMMM ())
BAM (1
RMMM 80 (3))
SA1 (-3)
MEE (ME)
NAA (PU)
MART
MM (11)
4N40-1
(-1)
4-1)
RRMMI te
4
P5A2 W-4)
ilpor
14142
(63)
C RRMM21 (RMM2
CANO (15)
800
BAR
AA @PMA
(t)
BR MME (U)
SAAI
0SAM 15
@R BRMM2
SA2 (1)
Figure 3: Results on JR55. Up: HiBaNG discovered causal links between indices on JR55 where the weights are from the Bernoulli posterior of the graph (links with weights less than 0.2 are not shown) and thicker links indicate stronger connections. Down: PCMCI+ discovered causal links.
5.3 Qualitative Analysis on Climate Reanalysis Data
We qualitatively analyze our method's performance on climate data obtained from the Japanese reanalysis of the atmosphere (JRA55) (Kobayashi et al., 2015), detailed in Soction C We compare HiBaNG and PCMCI+ with lag Tmax 6 in Figure We qualitatively explain the physical interpretability of relationships discovered by PCMCI and HiBaNG: 1) ENSO (MET) autocorrelation HiBaNG captures autocorrelations in MEI extending back to & 5, aligning with the known persistence of ENSO over approximately 6 months, as discussed in Harries & O'Kane (2021) (Figure 5). PCMCI only captures MEI autocorrelations up to 1-3, suggesting more limited temporal sensitivity. 2) MEI->PS1 directionality: HiBaNG correctly infers MEI(t 1) P))P)())),) consistent with the established ENSO-to-PSAI influence via modulation of midlatitude tropospheric flow (O'Kane Franzke, 2025). PCMCI infers the reverse direction (PSA1 1)-+)II( 3)
PNA links: Both methods identify the PNA autocorrelation at a one-month lag, consistent with observed behavior. However, HiBaNG additionally identifies the MEI(t-1) (EI((t)))-PNA(t edge, also seen in Harries & O'Kane (2021) (Figure 9b), supporting known teleconnections. PCMCI instead finds NAO (1 1) (L)((,, which is plausible but not as directly supported by prior work. 4) Northern Hemisphere modes: HiBaNG provides a more complete representation of the interconnected dynamics among AO, NAO+, NAO-, and AR. For example, the HiBaNG graph includes: autocorrelations at 1-1. These are consistent with Harries & O'Kane (2021) (Figures 6 and 7), and are not captured by PCMCI. Overall, the HiBaNG graph shows larger and more diverse set of edges consistent with those previously inferred from the JRA55 reanalysis via Bayesian structure learning as reported in qualitative demonstration
-1N((t-))NNN+++)) PSA(())+AR();; and NAO+ and NAO-
Harries & O'Kane )2021). We note that this experiment is intended primarily of our framework's 's interpretability and practical applicability, with PCMCT serving as a reference. We do not claim that HiBaNN should systematically outperform PCMCI on this task, particularly since PCMCI is more expressive in principle due to its support for nonlinear conditional independence tests.
6 Conclusion
We have presented novel Bayesian VAR model tailored to Granger causal discovery on MTS data in low-data settings. Our method leverages hierarchical Bayesian framework that separates Granger causal relationships into binary causal graphs and real-valued weights. Through extensive experiments on synthetic, semi-synthotic, and real-world datasets, we have demonstrated that our approach can perform better in low-data regimes. Regarding limitations, our method based on the Granger causality framework, which assumes that causal relationships are reflected in timo-laggod dependencics. While this assumption appropriate for many applications, it may not hold universally, and practitioners should exercise caution when interpreting the results.
12

---

Published in Transactions on Machine Learning Research (02/2026)
Broader Impact Statement
This paper introduces a novel method for discovering causal relationships from observational data, grounded in the framework of Granger causality and a factorized representation of causal structure. By leveraging these assumptions, the method enables scalable and interpretable causal discovery, which can benefit applications in fields such as oconomics, neuroscience, and climate science, where temporal data is abundant. The approach opens up new possibilities for causal inference in high-dimensional settings, providing a loundation for future work that can relax or adapt these assumptions to broader domains.
Because the method relies assumptions such as Granger causality and factorized causal structures, there is a risk of drawing incorrect conclusions if these assumplioss do not hold. Misinterpretation misuse of inferred causal relationships could lead to flawed decisions or reinforce biases present in the data, especially in high-stakes domains.
Acknowledgments
The authors thank the reviewers and action editor for the insightful reviews.
References
Daniel Felix Ahelegbey, Monica Billio, and Roberto Casarin. Sparse graphical vector autoregression: A Bayesian approach. Arnals of Economics and Statislics/ Statislics/Amaale Annales d' d'Economie et de Statislique, (123/124), 333 333-361, 2016.
David Alvarez Melis and Tommi Jaakkola. Towards robust interpretability with self-explaining neural networks. Advances in neural information processing systems, 31, 2018.
Yashas Annadani, Nick Pawlowski, Joel Jennings, Stefan Bauer, Cheng Zhang, and Wenbo Gong, BayesDAG: Gradient-basse posterior inference for causal discovery. Advances in Neural Information Processing Systems, 36, 2023.
Charles K Assaad, Emilie Devijver, and Eric Gaussier. Survey and evaluation of causal discovery methods for time series. Journal of Artificial Intelligence Research, 73:767-819, 2022.
Nicolas Bacaër and Nicolas Bacaër. Lotka, Volterra and the predator prey system (1920- 1926). short history of mathematical population dynamics, pp. 71-76, 2011.
Yoshua Bengio, Salem Lahlou, Tristan Deleu, Edward Hu, Mo Tiwari, and Emmanuel Bengio. GFlowNet foundations. Journal of Machine Learning Research, 24(210):1 55, 2023.
Yoav Benjamini and Yosef Hochberg. Controlling the false discovery rale: A practical and powerful approach to multiple testing. Journal of the Royal statistical society: series B (Methodological), (7(1):289-300 1995.
Monica Billio, Roberto Casarin, and Luca Rossini. Bayesian nonparametric sparse models. Journal of Econometries, 212 2019.
Konstantina Biza, Ioannis Tsamardinos, and Sofia Triantafillou. Tuning causal discovery algorithms. In International Conference on Probabilistic Graphical Models, pp. 17-28. PMLR, 2020.
Edwin V Bonilla, Pantelis Elinas, He Zhao, Maurizio Filippone, Vassili Kitsios, and Terry O'Kane, Variational dag estimation via state augmentation with stochastic permutations. arXio preprint arXiv:2402.0664 02644, 2024.
Jörg Breitung and Norman R Swanson. Temporal aggregation and spurious instantaneous causality in multiple time series models. Journal of Time Series Analysis, 23(6)6651-665, 2002.
Bart Bussmann, Jannes Nys, and Steven Latré Neural additive vector autoregression models for causal discovery in time series. In International Conference on Discovery Science, pp. 446-460. Springer, 2021.
13

---

Published in Transactions on Machine Learning Research (02/2026)
John Canny. GaP: factor model for discrete data. In Proceedings of the 27th annual international ACM SIGIR conference on Research and development in information retrieval, pp. 122-129) 2004.
George Casella and Edward I George. Explaining the gibbs sampler. The Americon Statistician, 46(3): 167-174, 1992.
Yuxiao Cheng, Runzhao Yang, Tingxiong Xiao, Zongren Li, Jinli Suo, Kunlun He, and Qionghai Dai. CUTS: Neural causal discovery from irregular time-series data. In The Eleventh International Conference on Learning Representations, 2023.
Yuxiao Cheng, Lianglong Li, Tingxiong Xiao, Zongren Li, Jinli Suo, Kunlun He, and Qionghai Dai. CUTS+: High-dimensional causal discovery from irregular time-series. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 38, pp. 11525-11533, 2024a.
Yuxiao Cheng, Ziqian Wang, Tingxiong Xiao, Qin Zhong, Jinli Suo, and Kunlun He. CausalTime: Realistically generated time-scrios for benchmarking of causal discovery. In The Twelfth International Conference on Learning Represenlations, 2024b.
Chris Cundy, Aditya Grover, and Stefano Ermon. BCD Nets: Scalable variational approaches for Bayesian causal discovery. Advances in Neural Information Processing Systems, 34:7095-7110, 2021.
Thomas Dean and Keiji Kanazawa. A model for reasoning about persistence and causation. Computational intelligence, 5(2):142 150, 1989.
Tristan Deleu, Mizu Nishikawa- Toomkyy Jithendaraa Subramanian, Nikolay Malkin, Laurent Charlin, and Yoshua Bengio. Joint. Bayesian inference graphical structure and parameters with a single generative flow network. Advances in Neural Information Processing Systems, 36, 2023.
Selva Demiralp and Kevin D Hoover. Searching for the causal structure of vector autoregression. Orford Bulletin of Economics and statistics, 65:745-767, 2003.
Chenchen Fan, Yixin Wang, Yahong Zhang, and Wenli Ouyang. Interpretable multi-scale neural network for Granger causality discovery. In ICASSP 2023-200 2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 1-5. IEFE, 2023.
Thomas Ferguson. Bayesian analysis of some nonparametric problems. The anmals of statistics, pp- 209 230, 1973.
Emily Fox, Erik B Sudderth, Michael I Jordan, and Alan S Willsky, Bayesian nonparametric inference of switching dynamic lincar models. JEEE Transactions on signal processing, 59(4):1569-1585 2011.
Karl J Friston, Lee Harrison, and Will Penny. Dynamic causal modelling. Neuroimage, 19(4):1273 1302, 2003.
Tomas Geffner, Javier Antoran, Adam Foster, Wenbo Gong, Chao Ma, Emre Kiciman, Amit Sharma, Angus Lamb, Martin Kukla, Nick Pawlowski, et al. Deep end-to-end causal inference. arXio preprinl arXiv:2202.02111 2022.
Edward I George, Dongchu Sun, and Shawn Ni. Bayesian stochastic search for var model restrictions. Journal of Econometrics, 142(1):5534 (155))558880 580, 2008.
Charles J. Geyer. The K-truncated Poisson distribution. URL https://www.sst umn. edu/geyer/aster Librarry/stter/bbc/kkp... pdf]
Satyajit Ghosh, Kshitij Khare, and George Michailidis. High-dimensional posterior consistency in Bayesian vector autoregressive models. Journal of the American Statistical Association, 2018.
Satyajit Ghosh, Kshitij Khare, and George Michailidis. Strong selection consistency of Bayesian vector autoregressive models based on a pseudo-likelihon approach. The Annals of Statistics, 49(3) 2021.
14

---

Published in Transactions on Machine Learning Research (02/2026)
Michael B Giles. Algorithm 955: Approximation of the inverse Poisson cummlative distribution function. ACM transactions mathemalical softwwore (TOMS), 42(1))1-22, 2016.
Walter R Gilks, Sylvia Richardson, and David Spiegelhalter. Markow chain Monte Carlo in practice. CRC press, 1995.
Clark Glymour, Kun Zhang, and Peter Spirtes. Review of causal discovery metuangods based on graphical models. Frontiers in genetics, 10:524, 2019.
Chang Gong, Di Yao, Chuzhe Zhang, Wenbin Li, Jingping Bi, Lun Du, and Jin Wang. Causal discovery from temporal data. In Proceedings of the 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, pp. 5803 5804, 2023.
Wenbo Gong, Joel Jennings, Cheng Zhang, and Nick Pawlowski. Rhino: Deep causal temporal relationship learning with history-deeenddet noise. arXiv preprint urXiv:2210.14706, 2022.
Prem K Gopalan, Laurent Charlin, and David Blei. Content-based recommendations with Poisson factorization. Advances in neural information processing systems, 27, 2014.
Clive WJ Granger. Investigating causal relations by econometric models and cross-spectral methods. Economelrica: Journal of the Econometric Sociely, pp.
1969.
Chuan Guo, Gcoff Pleiss, Yu Sun, and Kilian Q Weinberger. On calibration of modern neural networks. In International conference on machine learning, pp. 1321-1330. PMLR, 2017.
Dylan Harries and Terence J O'Kane. Dynamic Bayesian networks for evaluation of Granger causal rela- tionships in climate reanalyses. Journal of Advances in Modeling Earth Systems, 20020M540000 2021.
Biwei Huang, Kun Zhang, Jiji Zhang, Joseph D Ramney, Ruben Sanchez-Romes Clark Glymour, and Bernhard Schölkopf. Causal discovery from neterogencouu//nonstatiinn data. J. Mach. Learn. Res., 21 (89):1 53, 2020.
Aapo Hyvärinen, Kun Zhang, Shohei Shimizu, and Patrik O Hoyer. Estimation of a structural vector autoregression model using non-gaussianity. Journal of Machine Learning Reseorch, 11(5). 2010.
Guido W Imbens and Donald Rubin. Causal inference in statistics, social, and biomedical sciences. Cambridge university press, 2015.
Eric Jang, Shixiang Gu, and Ben Poole. Categorical reparamoterization with Gumbel-Softmax. arXiv preprint arXiv::611. 161100141, 2016.
Rahi Kalantari and Mingyuan Zhou. Graph gammn process linear dynamical systemss. In International Conference on Artificial Intelligence and Statistics, pp. 4060-4068. PMLR, 2021.
Rahi Kalantari, Joydeep Ghosh, and Mingyuan Zhou. Nomparametric Bayesian sparse graph linear dynamical systems. In International Conference on Artificial Intelligence and Statistics, pp. 1952-1960 PMLR, 2018.
Alireza Karimi and Mark R. Paul. Extensive chaos in the lorenz-96 model. Chaos: An interdisciplinary journal of nonlinear science, 20(4), 2010.
Saurabh Khanna and Vincent YF Tan. Economy statistical recurrent units for inferring nonlinear Granger causality. In International Conference Learning Representations, 2020.
Shinya Kobayashi, Yukinari Ota, Yayoi Harada, Ayataka Ebita, Masami Moriya, Hirokatsu Onoda, Kazutoshi Onogi, Hirotaka Kamahori, Chiaki Kobayashi, Hirokazu Endo, et al. The jra-55 reanalysis: General specifications and basic characteristics. Journal of the Meteorological Society of Japan Ser. II, 93(1):5-48, 2015.
Ananya Kumar, Percy S Liang, and Tengyu Ma. Verified uncertainty calibration. Advances in Neural Information Processing Systems, 32, 2019.
15

---

Published in Transactions on Machine Learning Research (02/2026)
Robert B Litterman. Forecasting with Bayesian vector autoregressions five years of experience. Journal of Business & Economic Slatistics, 4(1):25 38, 1986.
Jeremiah Liu, Zi Lin, Shreyas Padhy, Dustin Tran, Tania Bedrax Weiss, and Balaji Takshminarayanan. Simple and principled uncertainty estimation with deterministic deep learning via distance awareness. Advances in neural information processing systems, 33:7498-777 7512, 2020.
Lars Lorch, Jonas Rothfuss, Bernhard Schölkopf, and Andreas Krause. DiBS: Differentiable Bayesian structure learning. Advances in Neural Information. Processing Systems, 34:24111-244 1-1122,,,, 2021.
Edward N Lorenz. Predictability: A problem partly solved. In Proc. Seminar on prediclability, volume 1. Reading, 1996.
Sindy Löwe, David Madras, Richard Zemel, and Max Welling. Amortized causal discovery: Learning to infer causal graphs from time-scries data. In Conference on Causal Learning and Reasoning, pp. 509-525. PMLR, 2022.
Aurélie C Lozano, Naoki Abe, Yan Liu, and Saharon Rosset. Grouped graphical Granger modeling for gene expression regulatory networks discovery. Bioinformatics, 25(12):111--1118, 2009.
Helmut Lütkepohl. New introduction to multiple time series analysis. Springer Science & Business Media, 2005.
Damian Machlanski, Spyridon Samothrakis, and Paul Clarke. Hyperparameter tuning and model evaluation in causal effect estimation. arXiv preprint arXiv:2303.014.3 2023.
Chris J Maddison, Andriy Mnih, and Yee Whye Teh. The concrete distribution: A contimuous relaxation of discrete random variables. arXiv preprint arXi::6611.00712 2016.
Ričards Marcinkevičs and Julia E Vogt. Interpretable models for Granger causality using self-cxplaining neural networks. In International Conference on Learning Representations, 2021.
V Mihajlovic and Milan Petkovic. Dynamic Bayesian networks: A state of the art Universily of Twenle Document Repository, 2001.
Silvia Miranda-Agrippinn and Giovanni Ricco. Bayesian vector autoregressions: Estimation. In Oxford Research Encyclopedia of Economics and Finance. 2019.
Alessandro Montalto, Sebastiano Stramaglia, Luca Facs, Giovanni Tessitore, Roberto Prevete, and Danicle Marinazzo. Neural networks with non-uniform embedding and explicit validation phase to assess Granger causality. Neural networks, 71:159-174 1711, 2015.
Kevin Patrick Murphy. Dynamic Bayesion networks: Representation, inference and learning. University of California, Berkeley, 2002.
Jouchi Nakajima and Mike West. Bayesian analysis of latent threshold dynamic models. Journal of Business & Economic Statistics, 31(2):151-164, 2013.
Meike Nauta, Doina Bucur, and Christin Seifert. Causal discovery with attention-based convolutional neural networks. Machine Learning and Knowledge Extruction, 111):19, 2019.
William B Nicholson, David S Matteson, and Jacob Bien. VARX-L: Structured regularization for large vector autoregressions with exogenous variables. International Journal of Forecasting, 33(3):627 651, 2017.
NumPyro. Truncated and folded distributions. truncated distributions... bxuucateeddisstriutttiosss. html
URL https://num.pyett http:///w.ppprro.a///////ttttttt
Junier Oliva, Barnabás Póczos, and Jeff Schneider. The statistical recurrent unit. In International Conference Machine Learning, pp. 2671-17 -2600. PMLR, 2017.
16

---

Published in Transactions on Machine Learning Research (02/2026)
Terence J. O'Kane and Christian L.E. Franzke. Chapter 8 Pacific-South American pattern. In Bin Guan (ed.), Almospheric Oscillations, pp. 171-181 Elsevier, 2025. ISBN 978-0-443-15638-00 doi: httpss//doi.org/ 10.1016/13978-0-10 13978-0--7------ 16638-0.00008-3 URL http://wwwceeccccccecccc////ss/s
B97804441163300000333
Roxana Pamfil, Nisara Sriwattanaworachai, Shaan Desai, Philip Pilgerstorfer, Konstantinos Goorgatzis, Paul Beaumont, and Bryon Aragam. DYNOTEARS: Structure learning from time-series data. In International Conference on Artificial Intelligence and Statistics, pp. 1595-1605 PMLR, 2020.
Judea Pearl. Causal inference in statistics: An overview. Stalistics surveys, 3:96 146, 2009.
Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Causal inference on time series using restricted structural equation models. Advances in Neural Information Processing Systems, 26, 2013.
Jakob Runge. Causal network reconstruction from time series: From theoretical assumptions to practical estimation. Choos: An Interdisciplinary Journal of Nonlinear Science, 28(7):075310, 2018.
Jakob Runge. Discovering contemporaneons and lagged causal relations in autocorrelated nonlinear time series datasets. In Conference on Uncertainty in Artificial Intelligence, pp. 1388-1397. PMTRR, 2020.
Jakob Runge, Peer Nowack, Marlene Kretschmer, Seth Flaxman, and Dino Sejdinovic. Delecting and quantifying causal associations in large nonlinear time series datasets. Science advances, 5(11) ((11))eaau4996, 2019.
Aaron Schein, Scott Linderman, Mingyuan Zhou, David Blei, and Hanna Wallach. Poisoon-andommzzdd gamma dynamical systems. Advances in Neural Information Processing Systems, 32, 2019.
S Seabold and PJ Statsmodels. Econometric and statistical modeling with python. In Proceedings of the Ith Python in science conference, pp. 57-61, 2010.
Pedro Shiguihara, Alneu De Andrade Lopes, and David Mauricio. Dynamic Bayesian network modeling. learning, and inference: A survey. IEEE Access, 9:117639 117648, 2021.
Ali Shojaie and Emily B Fox. Granger causality: A review and recent advances. Annual Review of Statistics and Its Application, 9:289-30 319, 2022.
Stephen Smith, Karla Miller, Gholamreza Balimi-Khorshin Matthew Webster, Christian F Beckmann, Thomas E Nichols, Joseph Ramsey, and Mark W Woolrich. Network modelling methods for fmri. Neuroimage, 54(2):875-891, 2011.
Peter Spirtes, Clark Glymour, Richard Scheines, and David Heckerman. Causalion, prediction, and search. MIT press, 2000.
Norman R Swanson and Clive WJ Granger. Impulse response functions based on a causal approach to residual orthogonalization vector autoregressions. Journal of the American Statistical Association, 92 (437))357-367, 1997.
Alex Tank, Ian Covert, Nicholas Foti, Ali Shojaic, and Emily Fox. Neural Granger causality for nonlincar time series. stat, 1050:16, 2018.
Ryan Thompson, Edwin V Bonilla, and Robert Kohn. Contextual directed acyelic graphs. In International Conference on Arlificial Intelligence and Statislics, pp. 2872 2880. PMLR, 2024.
Ryan Thompson, Edwin V Bonilla, and Robert Kohn. Prodag: Projected variational inference for directed acyelic graphs. In The Thith Annual Conference on Neural Information Processing Systems, 2025.
Alexander Tong, Lazar Atanackovic, Jason Hartford, and Yoshua Bengio. Bayesian dynamic causal discovery. In A causal view dynamical systems, NewrIPS 2022 workshop, 2022.
Sumanth Varambally, Yian Ma, and Rose Yu. Discovering mixtures of structural causal models from time series data. Tn International Conference Machine Learning, 2024.
17

---

Published in Transactions on Machine Learning Research (02/2026)
Yueming Wang, Kang Lin, Yu Qi, Qi Lian, Shaozhe Feng, Zhaohui Wu, and Gang Pan. Estimating brain connectivity with varying-length time lags using recurrent neural network. IEEE Transactions on Biomedical Engineering, 65(9):1953 1963, 2018.
Robert L Wolpert, Merlise A Clyde, Chong Tu, et al. Stochastic expansions using continuous dictionaries: Lévy adaptive regression kernels. The Annals of Statisties, 39(4):196-1196, 2011.
Tomasz Wožniak. Bayesian vector autoregressions. Australion Economic Review, 49(3):365-23) 2016.
Dongxia Wu, Tsuyoshi Idé, Georgios Kollias, Jiri Navratil, Aurelie Lozano, Naoki Abe, Yian Ma, and Rose Yu. Learning Granger causality from instance-wise self-attentive Hawkes processes. In International Conference on Artificial Intelligence and Statistics, pp. 415-423. PMLR, 2024.
Tailin Wu, Thomas Breuel, Michael Skuhersky, and Jan Kautz. Discovering nonlinear relations with minimum predictive information regularization. arXiv preprint arXiv:2001.01883 2020.
Jaewon Yang and Jure Leskovec. Community-afllittionn graph model for overlapping network community detection. In 2012 TEFE 12th international conference on data mining, pp- 1170-1174 IEFE, 2012.
Jacwon Yang and Jure Leskovec. Structure and overlaps of ground-truth communities in networks. ACM Transactions Intelligent Systems and Technology (TIST), 5(2):1 35, 2014.
Ming Yuan and Yi Lin. Model selection and estimation in regression with grouped variables. Journal of the Royal Statistical Society Series B: Statistical Methodology, (8(1):49-67, 2006.
He Zhao, Lan Du, and Wray Buntine. Leveraging node attributes for incomplete relational data. International conference on machine learning, pp. 4072-4081. PMLR, 2017.
He Zhao, Lan Du, Wray Buntine, and Mingyuan Zhou. Dirichlet belief networks for topic structure learning. Advances in neural information processing systems, 31, 2018a.
He Zhao, Lan Du, Wray Buntine, and Mingyuan Zhou. Inter and intra topic structure learning with word embeddings. In International Conference on Machine Learning, pp. 5892 5901. PMLR, 2018b.
MingYuan Zhou. Infinite edge partition models for overlapping community detection and link prediction. In AISTATS, pp. 1135-1143, 2015.
Mingyuan Zhou and Lawrence Carin. Negative binomial process count and mixture modeling. IEEE Transactions on Pattern Analysis and Machine Intelligence, (7(2):007-320, 2013.
Mingyuan Zhou, Lauren Hannah, David Dunson, and Lawrence Carin. Beta-negative binomial process and Poisson factor analysis. In Artificial Intelligence and Statistics, pp. 1462 1471. PMLR, 2012.
Wanqi Zhou, Shuanghao Bai, Shujian Yu, Qibin Zhao, and Badong Chen. Jacobian regularizer-based neural Granger causality, In International Conference on Machine Learning, 2024.
A Inference via Gibbs Sampling
Sampling A With the conjugacy of normal distributions, one can sample the entries of AT one by one by:
(O,,
if G =0
AJ,
(7)
(NT,
otherwise
18

---

Published in Transactions on Machine Learning Research (02/2026)
where:
T
o =
[x.GT. Σ
2
+
Wils
(8)
tml
Σ Fil
(9)
Tiaax
Fist
Σ 1'
G.
G.
(10)
Sampling Vi, and A With the conjugacy between normal and gamma distributions, one can sample क(,,
and from their conditional gamma posteriors:
ViJ
Gamma (1.5, ( 1))
(11)
-1
Tinax
N
A Gamma
1+T/2,
1+
Σ
A G(
(12)
j=1
Sampling M
With GBPL, we can sample:
(
1)
Categoricaly
Th )ur (..)
if G =0
Mr
(13)
DV 2'
TPoissony (9T,).
otherwise
where di,- Ek To
Sampling M.
With the relationships between Poisson and multinomial distributions, we can sample:
=== LT.j. ..... ..... 1).
Fiju
[..., I..,m......... Multinomialy mī.ji
(14)
where J. OT,NTO
Sampling H T7 and rg
N
0 Gamma
of +
Σ
MT
(15)
Light
di. + r Σラ 03,4
j=1
N
④ Gamma
6j
Mij.k
MT
(16)
of + r Σm 07
1
rk Gamma 1/K+
Σ
Mij
(17)
с + Li=1 03.0
N
1=1,==
19

---

Published in Transactions on Machine Learning Research (02/2026)
Table 3: Hyperparameter settings.
Tiocs
# hidden layers
units
# training opochs
Learning rate
Parameter space
VAR (FBH)
(1,3,5)
NA
NA
NA
NA
NA
NA
BVAR(d)
(1.3.5)
NA
NA
NA
NA
NA
NA
BVAR(e)
(13355
NA
NA
NA
NA
NA
[0.011 0.05
SRU
NA
10
2000
1.0o-3
M2
[0.01, 1.0)
μ4 (0.01,
<SRU
NA
2
10
2000
1.0e-3
50
(0.01, 0.00
Pa [0.01, 1.0]
[0.0, 3.00
GVAR
(1,3,5)
2
50
1,000
1.0c-4
7 [0.0,
HiRaNG
(1.3.5
NA
NA
10,000
NA
NA
(1,2,3)
Sampling of and 67 By introducing auxiliary variables from the Chinese Restaurant Table (CRT)
distribution Zhou et al. 2012 Zhou & Carin
2013
sample:
1 CRT
Miskiaf
(18)
d Gamma 1+ 1+
(19)
1+ K
log(1 + r E5 /
k=1
k-
N
oJN ~ CRT
M bj
(20)
K 83 Gamma
1+ 1+Z0v
1
(21)
1 +
K Σ
log(1 + rg ΣM /
k=1
Sampling dy of and e
1
di ~ Gamma
Σ
0 + 1
(22)
Σ
0%
+
N
1
6 Gamma
6 +
1,
(23)
Σ
+1
1
C Gamma
2,
(24)
r +1
B
Synthetic Data
We conduct experiments on synthetic data generated from a. VAR model specified in Eq. 1 to test whether our method can discover the ground-truth graphs. Given N 16 and Tmax 6, construct {AT) Tam by first specifying the nonzero entries (i.e., the ground-truth causal graphs) and then, for each nonzero entry, sampling A[ (0.10.) We then generate T = 1,000 samples, initialize from standard normal distribution, and sample ~ N (0, 0.01). We show the results of our method in Figure 4 where we also fit. a randomly initialized VAR (Seabold Statsmodels] 2010) to the data as a reference. The ground-trute -tunh graphs different lags exhibit diverse patterns, and our method recovers them well (also reflected by higher AUROC and AUPRC). Unlike VAR, our method directly infers binary graphs without using thresholds hypothesis tests.
20

---

Published in Transactions on Machine Learning Research (02/2026)
Algorithm 1: Inference Algorithm for HiBaNG
input :MTS data X, number of lags Tmax, hyperparameter V output: Posterior samples of {A'] ymma and {G") Tmax Initialize all the variables; while Not converyed do for i= 1...... . do
| Sample A;; end
for
= 1...... Tmaz do for i 1... N. j...... NN do Sample M and Mij.i end Sample c²;
1 do
| Sample af and ;;; end
for l= K do Sample di c Th:
end
for
i= I. N do
for k = 1=...KK K do Sample o and of
end end for
== 1= do
for
1...NN N do Sample Aigs Vhy Gī
end
end
end
end
C Further Introduction to the Climate Reanalysis Data
Following Harries & O'Kane 2021
compute 13 indices that diagnose the activity of the major atmospheric, tropospheric, and convective global climate modes at monthly resolution from 1960 to 2005, resulting in an MTS dataset with N = 13 and T = 551. The climate indices' names are shown in Table 7 Among SouthSn Oscillation Index (MEI) representative in that its time
the indices, the Multivariate series associated with regionally distributed, coherent responses the atmosphere and surface ocean. In Figure 5 (a,, we plot the MEI time series an example; it characterizes the El Niño/La Niña cycle. Positive values of the MEI are associated with El Niño periods, whereas negative values are associated with La Niña periods; the magnitude of the index is proportional to the strength of the event. For example, according to the MEI, there was a strong El Niño in April 1998. Figure 5]b) illustrates how much warmer (red) and cooler (blue) the surface air temperature was in that month with respect to the average April. This map shows a large warm patch over the eastern Pacific Ocean, which typical indicates that August 1988 fcatured a strong La Niña; the associated surface air temperature map shown in Figure 5(c) c). This map illustrates how different the surface air temperature was with respect to an average August. As is typical of La Niña events, the eastern Pacific region is anomalously cool.
Niño. Conversely, the MEI
21

---

Published in Transactions on Machine Learning Research (02/2026)
Table 4: Performance HiBaNG on FMRI different parameters.
Tmax
K
AUROC
AUPRC
1
1
5
0.711 + 0.044
0.550 + 0.037
1
1
20
0.725 + 0.046
0.552 0.034
1
1
50
0.714 + 0.041
0.550 + 0.035
1
2
5
0.715 + 0.051
0.544 0.037
1
2
20
0.713 + 0.052
0.543 0.036
1
2
50
0.712 + 0.048
0.543 0.033
1
5
0.705 + 0.058
0.539 0.035
1
20
0.710 + 0.060
0.539 0.038
1
50
0.706 + 0.044
0.538 0.034
3
5
0.737 + 0.039
0.590 0.038
3
20
0.731 + 0.032
0.585 0.037
3
50
0.739 + 0.036
0.591 0.038
3
5
0.732 + 0.031
0.583
0.039
3
20
0.728 + 0.012
0.597 0.014
3
50
0.724 + 0.036
0.584 0.044
3
3
5
0.723 + 0.028
0.578 0.040
3
3
20
0.731 + 0.034
0.582 = 0.042
3
3
50
0.736 + 0.034
0.585 = 0.044
5
1
5
0.736 + 0.022
0.589 = 0.023
5
1
20
0.742 + 0.029
0.594 = 0.026
5
1
50
0.742 + 0.030
0.593 0.027
5
2
5
0.732 4: 0.026
0.587 + 0.033
5
2
20
0.743 +: 0.030
0.592 + 0.031
5
2
50
0.740 +: 0.029
0.591 + 0.029
5
3
5
0.733 + 0.017
0.583 + 0.031
5
3
20
0.738 + 0.030
0.590 + 0.032
5
3
50
0.734 +: 0.023
0.589 + 0.035
D Empirical Computational Performance
In addition to the complexity analysis Section 3.3 we empirically study computational efficiency by examining two aspects: the running time of single Gibbs sampling iteration and the number of iterations required for convergence. For all experiments in this section,
method an Apple laptop with vary the
MI Pro processor. We report the running time (seconds per Gibbs sampling iteration) number of variables and time steps on Lorenz 96 in Figure 6c Inference is efficient and scales gracefully with T. Figure shows HiBaNG's MSE over iterations on Lorenz 96. Although we set the maximum number of Gibbs sampling iterations to 10,000, in most cases our method converges in around 200 iterations. We also observe that low-data regimes (c.g., T 100 on Lorenz 96), samples from our model may exhibit higher variance early in inference, as expected with fewer time steps. general, our method converges within half hour laptop in most cases.
22

---

Published in Transactions on Machine Learning Research (02/2026)
Table 5: SHD for VAR (FBH), GVAR, and HiBaNG with Tmax computed over replicates on each dataset.
Mcans and standard deviations are
L96
LV
FMRI
T' 100
500
T 200
= 500
VAR (FBH)
72.00
125.60 2-1.88
451.40++69.88
26.00=1.10
GVAR
373.20++39.77
93.60+555.23
147.80+ 80.36
111.00++39.78
51.00+16.79
HiBaNG
112.91
71.43+3.52
186. 16.113+21.88
49.73++4.72
25.17=1.43
Table 6: SHD for VAR (FBH), GVAR, and HiBaNG with Tmax 3. Means and standard deviations are computed over replicates cach dataset.
T.66
LV
FMRT
T 100
T = 500
T 200
= 500
VAR (FBH)
79.60#3.32
84.40:4..18
75.40
25.80+1..33
GVAR
577.40+47.51
133. 133.20:107.22
213.80::117.9
128.60
50.00+18.99
HiBaNG
114.75=1.68
73.1842.67
75.30+9.21
45.50
24.78
Table 7: Climate index names.
Index name
Abbreviation
Atlantic Oscillation
AO
Indian Ocean Dipole
IOD
Multivariate El Niño Southern Oscillation Index
MEI
North Atlantic Oscillation (positive and negative phases)
NAO+/-
Atlantic Ridge patterns
AR
Scandinavian blocking patterns
SCAND
Pacific North American patterns
PNA
Pacific South American patterns
PSA1/2
Southern Annular Mode
SAM
Wheeler-Hendon Madden-Julian oscillation
RMM1/2 1/2
23

---

Published in Transactions on Machine Learning Research (02/2026)
kamic
Prde aieroc 18, Sanic 1it
Fndd sbe, wwr.nc. 1-8, sami
Fedobe sterne 1/8 samic
Pauds sicacc Lot Dimic
strrnc 10, Dami ais
Fedeiiis dicani iit Jami
Figure 4: Synthetic dataset. First to third rows: ground-truth graphs, coefficients VAR, and GC graphs of ours at different lags. Columns: T T= =l 6. The AUCROC and AUPRC scores each lag of VAR and ours are shown in the sub-captions. Mean AUCROC over all the lags is 0.99 (VAR) and 1.0 (ours); Mean AUPRC is 0.93 (VAR) and 0.96 (ours).
m
nhn
mar nw Limp
ma In M Mr
yy
-2
-3
Niño
Niña
1960
1965
1970
1975
1980
1985
1990
1995
2000
2005
(a)
32TM
301
00'5
23
60'5
122/E
300'
120W
6PW
-15
:23/E
1.0'
3227W
-15
(b)
(c)
(d)
Figure 5: (a) MEI index. (b) The geographical locations of the indices. (c) Fl Niño temperature anomaly
April 1998. (d) La Niña temperature anomaly August 1988.
E More Analysis on Bayesian nonparametrics
To further analyze the shrinkage mechanism of our model with respect to the truncation level K, for each K. which we denote by R E K. For 5, 5, we compute r5 = modian R). Since = is small, we use it as a reference point and define the Kmmar
run (Figure 2 we obtain the matrix of factor weights {r [rkkk=
number of effective latent factors & as # > rōjk=1,
Under this definition, we have lg 5+rmax by construction. For K E (2), we report La 6x-6 lg averaged 2
five with different random seeds. This quantity measures how much the number of effective latent.
24

---

Published in Transactions on Machine Learning Research (02/2026)
2.0
6-120
1.8
1.4
1.4
1.2
1.0
0.8
100
350
400
650
500
100
230
15.00
Iteration
Beration
(a) L96, T 100
(b) T.96, T 500
(c) Running times
Figure 6: (a, b) HiBaNG's (V 2) empirical converge Lorenz 96. For better visualisation, we show
MSE in the iterations between [a, b] where is the first iteration that MSE goes below 2.0 and b 500. (c) Running times in seconds per MCMC iteration on on Lorenz 96.
Table Difference in the number of effective latent factors, fx - Es, for K {20, 50) across datasets and lag orders. The reference value is computod from the case
-x-
L96, T
100
L96, T
500
IV, T
200
LV, T 500
FMRI
K = 20
K 50
K = 20
K = 50
K = 20
K = 50
K 20
K = 50
K 20
K = 50
Timxx =1
0.000
-0.667
-0.667
-0.333
-0.333
0.000
1.667
4.000
0.000
0.000
Thuxx =3
0.667
2.000
2.667
2.000
2.333
0.667
2.000
2.667
0.000
2.667
Tmax 55
2.000
3.000
1.333
4.667
3.000
1.000
2.000
2.667
-1.333
0.667
factors deviates from the reference case K = 5, thereby indicating whether the model activates substantially more factors when the truncation level is increased.
We emphasize that la is not intended to be a precise or exhanstive measure of the true number of active latent factors, but het a coarse and interpretable proxy for assessing posterior shrinkage behavior. Since the notion of an "effective" factor inherently ambiguous in hierarchical Bayesian models with continuous weights, any such metric is necessarily approximate. Our goal therefore not exact quantification, but provide simple diagnostic that reflects whether increasing K leads a systematic increase in the number of activated components. As shown in Table R most values of fx - ls remain small across all datasets, lag orders, and truncation levels. Even when increasing K from 5 to 50, the number of effective latent factors grows only marginally and, in some cases, even decreases. This indicates that the posterior selectively activates only a limited number of latent factors, while the remaining capacity remains unused.
F
Convergence Analysis
For a detailed convergence analysis of our algorithm, we set V = 2, Tmax = 3, and K = 20, and run 20 independent Gibbs sampling chains with different random seeds and initializations on the FMRI dataset. Consistent with our other experiments, each chain is run for 10,000 iterations, of which the first 5,000 are discarded as burn-in. From the remaining iterations, we retain one sample every 10 itcrations, resulting 500 samples per chain.
We first compute the posterior edge probabilities p (mean and standard deviation) for each pair of nodes at cach lag by pooling the 500 post-burn-in samples across all chains. We then select a representative subset of 30 triplets (i,j.) (where i and j are node indices and is the lag index) with the largest posterior means for visualization. To assess convergence and mixing, we report the Gelman Rulman diagnostic (R) and effective sample size (ESS) for both G and A for the selected node pairs. Table reports the convergence diagnostics. For both G and A, Å and ESS are close to 1, indicating reasonably good agreement across the 20 independent chains. Moreover, the effective sample sizes consistently high and close to the total
R and ESS:
25

---

Published in Transactions on Machine Learning Research (02/2026)
number of retained samples, suggesting strong mixing and low autocorrelation. These results demonstrate that our Gibbs sampler achieves stable and reliable posterior inference for all selected interactions.
Table 9: Analysis of and A in multiple chains on FMRI.
G
A
p
R
ESS
R
ESS
10
15
0.9001 = 0.29988
1.0013
6380.5
-0.28935 + 0.70166
0.99994
9000.7
9
8
0.3615 = 0.48046
1.0002
7447.9
-0.11128
2.7951
0.99959
9616.7
10
5
0.1468 = 0.35392
1.0011
8421.7
-0.014161
2.5222
1.0005
8795.9
12
11
0.118 0.32262
1.0001
8672.4
0.043694 = 2.9617
0.99935
9278.9
6
13
0.0924 + 0.2896
1.0009
9453
0.012286 = 2.6168
1
9663.2
2
10
0.0753 = 0.26389
1.0001
9352
0..18597 2.7518
1
9383.4
14
15
0.074 0.26178
0.99969
9083.2
0.029438 = 3.5335
1.0003
9653.4
3
10
0.0725 = 0.25933
0.9995
8889.1
0.059261 2.7133
0.9996
9507.9
1
11
0.0576 + 0.233
1.0011
9375.3
00.0049886 + 2.8795
0.99953
9554.4
13
15
0.0548 + 0.2276
1.0009
8844.8
0.070477 = 2.8283
0.99949
9184.5
1
8
0.0529 = 0.22385
1.0002
9159.9
0.026641 = 3.1418
0.99943
9614
6
7
0.0515 = 0.22103
0.99983
9802.8
-0.024542 + 4.5388
1.0002
9473.8
13
6
0.0504 = 0.21878
0.99947
9634.5
-0.043866 + 3.4608
0.99991
9584
11
1
0.0462 = 0.20993
0.99996
9065.2
-0.022 + 2.7709
0.99933
9097.7
6
14
0.0424 = 0.20151
0.99949
9337
0.019762 = 2.6693
0.99984
9238.7
7
3
0.042 = 0.2006
0.99935
9489.2
0.039643 = 3.1857
0.99974
8766.6
2
4
0.041 = 0.1983
0.99992
8932.8
0.015108 + 4.9719
0.99973
9161.6
10
4
0.0366 = 0.18779
1.0003
9267
0.047841 = 3.3277
0.99952
9003.4
10
8
0.0351 = 0.18404
0.99998
9416.1
0.013438 + 2.571
0.99941
9730.2
2
11
0.0338 = 0.18072
1.0005
9180.9
0.057287 = 3.4617
0.99993
9076.5
3
14
0.0323 + 0.1768
1
9177.5
0.029325 = 3.4226
0.99954
9608.5
7
14
2
0.0294 0.16893
0.99944
9598.2
0.019568 = 4.3586
1.0002
8831.7
8
10
0.027 + 06209
1.0006
9494.6
-0.02755 2.5024
0.99996
9472
8
6
0.0266
0.16092
0.99939
9360.1
00.012545 + 3.1591
0.99953
9658.2
10
12
0.0249 0.15583
1.0004
9446.7
0.033767 = 3.3515
1.0004
9664.7
7
1
0.0245 + 0.1546
0.99925
9573.5
0..00686
11.404
0.99977
9853.3
11
13
0.0234 0.15118
1.0006
9555.4
0.016013 = 3.8499
1.001
9381.4
5
4
0.0233 0.15086
0.99982
9487.6
00.34036 + 3.6502
0.99944
8998.4
15
1
0.0224 0.14799
1.0004
9047.4
0.038182 + 4.2819
0.99974
9340.3
14
13
0.0223 0.14766
1.0012
9922.7
-0.044632 + 3.1499
0.99968
9062.1
Reliability diagram: We evaluate the probabilistic calibration of posterior edge probabilities using reliability diagrams with respect to the ground-truth graph on the FMRI dataset. Due to the extreme sparsity of true odges, we construct bins with approximately equal numbers edges obtain stable empirical estimates. For each bin, we plot the mean predicted probability against the corresponding empirical edge frequency. Each reliability diagram is accompanied by histogram of predicted probabilities to contextualize the calibration
26

---

Published in Transactions on Machine Learning Research (02/2026)
Reliability diagram
Empirical frequency
0.8
0.6
0.4
0.2
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
Mean predicted probability
Prediction distribution
200
150
Count
50
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
Predicted probability
Figure 7: Reliability diagrams.
1015,
Figh probutio
12
0
bal 1-13 legert
jeta lnger
lagel 100000
22
Dgg 0000002
100
200
400
Sanplo
Surrplo inceoo
Figure 8: Edge probability traces. Colors indicate different chains.
curve. As shown in Figure 7 most predictions concentrate near zero, reflecting the sparsity of the underlying graph, while the calibration curve closely follows the diagonal in the low-probability regime, indicating reasonable calibration where most of the probability mass lies. This aligns with the CE results reported in Table
Edge probability traces: To visually assess convergence and mixing of the Gibbs sampler, we include multi-chain trace plots for selected edgo-specific variables. Because the latent variables are matrix-valued, it is not. feasible to display traces for all entries. Instead, visualize traces for representative subset of
27

---

Published in Transactions on Machine Learning Research (02/2026)
triplets (1. r), selected previously. In Figure 8 for each selected triplet (4.j. r), we plot its edge probability across the collected 500 samples for all chains on the same axis. These trace plots provide a qualitative diagnostic of convergence, illustrating that all chains explore the same stationary distribution, exhibit no visible nonstationarity, and mix adequately. We complement these visual checks with quantitative convergence diagnostics shown Table g
28


Reviews:
[{"review_id": "r5LnyO6RVQ", "reviewer": "Reviewer_gDsH", "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": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper introduces **HiBaNG**, a hierarchical Bayesian nonparametric framework for Granger causal discovery in multivariate time series, with a particular focus on low-data regimes where deep learning–based VAR models often struggle. The key contribution is an explicit probabilistic separation between binary Granger causal graphs and real-valued VAR coefficients, enabled through a Poisson-factorized hierarchical prior and a generalized Bernoulli–Poisson link function. This design allows structured sparsity, uncertainty-aware inference, and tractable Gibbs sampling. The methodology is carefully developed and theoretically motivated, and the experimental evaluation, especially on synthetic and semi-synthetic benchmarks, is thorough and convincing. At the same time, the approach relies on a stylized linear VAR assumption, which might not be useful in real-world applications. It introduces modeling choices (e.g., lag-wise independent graph priors and hierarchical Poisson constructions) whose necessity and generality would benefit from further clarification and justification.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The main claim of the paper, namely that the proposed hierarchical Bayesian nonparametric framework enables effective Granger causal discovery in low-data regimes, is supported by convincing empirical evidence. The authors present a comprehensive experimental evaluation on synthetic, semi-synthetic, and real-world datasets. In particular, the synthetic and semi-synthetic experiments are well controlled and directly test the stated advantages of the method in low-sample settings. That said, while the empirical results are strong, several modeling choices and claims would be further strengthened by more meticulous theoretical analysis and justification, particularly regarding temporal assumptions and prior construction.", "audience_yn": "Yes", "audience_explanation": "Yes, this paper is of high significance to the causal discovery, VAR, Granger causality, and time-series community. Granger causality is also tied to root cause analysis/diagnosis for systems with time-series data. Therefore, researchers working on anomaly detection would also benefit from this analysis. However, the restriction of the models to linear statistical methods with closed-form might limit their applicability to only controlled synthetic datasets. The authors need to work on the extension of this work to Deep methods for supporting real-world industrial applications.", "requested_changes": "**(1) Independence of Granger graphs across time lags**\n\nIn the multiple-lag setting, the paper states that “we have a separate generative process for the GC graph at each lag $\\tau$.” This appears to imply that the prior factorizes as\n$$\np(G^1, \\dots, G^{\\tau_{\\max}}) = \\prod_{\\tau=1}^{\\tau_{\\max}} p(G^\\tau),\n$$\ni.e., that the Granger causal graphs at different lags are *a priori* independent. Please clarify whether this interpretation is correct. If so, the authors should provide a clearer mathematical justification for why independent priors across time lags are reasonable, especially given that many real-world dynamical systems exhibit consistent or decaying causal effects across multiple delays. A short theoretical argument or discussion explaining why likelihood-based posterior coupling alone is sufficient (or preferable) would strengthen the modeling rationale.\n\n**(2) Linear VAR assumption vs. Deep VAR models**\n\nThe proposed framework is built on a linear Gaussian VAR likelihood, which is a stylized assumption that is rarely satisfied in complex real-world systems. In contrast, much of the recent state of the art in Granger causal discovery relies on nonlinear or neural-network-based VAR models (Deep VAR). While the paper provides an extensive and thoughtful literature review of Deep VAR methods, it remains unclear whether HiBaNG should be viewed primarily as (i) an interpretable simplification of Deep VAR approaches, or (ii) a fundamentally different modeling paradigm with distinct theoretical advantages.\n\nAt a minimum, the paper would benefit from a dedicated discussion or section outlining how the proposed hierarchical Bayesian graph prior could be extended to nonlinear or neural VAR settings, for example by replacing the linear predictor\n$$\nx_t = \\sum_{\\tau} (A^\\tau \\odot G^\\tau)x_{t-\\tau}\n$$\nwith a parameterized nonlinear function while retaining the same probabilistic treatment of $G^\\tau$. This section should include an appropriate mathematical formulation. Additionally, 1–2 sentences should be added to the Deep VAR portion of the related work to explicitly articulate how HiBaNG differs in modeling assumptions and inference goals.\n\nRelatedly, please clarify whether any preprocessing or transformations were applied to the real-world climate data to justify a linear VAR assumption. If no such preprocessing was used, does the paper implicitly claim that a linear VAR provides an adequate fit to the raw climate indices? This point requires clarification.\n\n**(3) Justification for the hierarchical Bayesian graph prior**\n\nThe motivation for adopting a hierarchical Bayesian construction over Granger causal graphs remains insufficiently justified. It is unclear what would fail, theoretically or empirically, if a simpler prior were used, such as\n$$\nG^\\tau_{i,j} \\sim \\mathrm{Bernoulli}(\\pi), \\quad \\text{or} \\quad \\pi \\sim \\mathrm{Beta}(a,b).\n$$\nThe authors should provide theoretical insight or formal arguments explaining why such Bernoulli (or Beta–Bernoulli) graph models are inadequate in this setting, and why the proposed Poisson-factorized hierarchical construction is necessary. Without this justification, the methodology risks appearing as a collection of sophisticated statistical components rather than a principled modeling choice grounded in theory.", "broader_impact_concerns": "No concerns about the ethical implications of this work. The authors have done a good job mentioned a detailed \"Broader Impact Statement\".", "additional_comments": "This is a well-written paper with clearly stated assumptions and thoughtful justifications supporting its main claims. The statistical methodology is carefully presented, and the experimental evaluation, particularly on synthetic and semi-synthetic datasets, is comprehensive and well executed. That said, the paper has several conceptual and modeling limitations, particularly regarding temporal assumptions, linearity, and the necessity of hierarchical graph priors, which are detailed in the requested changes above.", "extra_scores": {}, "extra_text": {}}, {"review_id": "wBNfoFKdaD", "reviewer": "Reviewer_Y2nF", "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": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "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 addresses the problem of discovering Granger causality among variable in multivariate times-series using Bayesian nonparametric framework. Focusing on the low-data regime, they propose a sparsity-encouraging prior for the underlying causal graphs and develops an Gibbs sampling algorithms for inference.\nThe proposed method is evaluated and  compared to several baseline for synthetic and real data, but there is not theoretical results.\n\nThe paper makes the following main claims:\n1) The proposed prior is claimed to enable \"interpretable, uncertainty aware inference\"\n2) The resulting algorithm is claimed to improve accuracy and calibration\n3) The resulting algorithm reduces overfitting\n4) The resulting algorithm has a reduced number of tunable hyperparameters.\n\nGenerally, the paper is well-written, well-structured, and easy to follow. However, section 3 and 5 could be significantly improved, see comments below.", "claims_and_evidence_yn": "No", "claims_explanation": "Ideally I would have answered: \"Yes to no claims\" and \"no for other\" above. I'll elaborate in the following:\n\nRegarding claim 1)\nThe abstract, introduction, and conclusion highlight interpretability in various ways, but it is completely ignored in the rest of the paper. Therefore, this claim is not supported. Moreover, it is somewhat unclear what \"interpretable, uncertainty-aware inference\" even means.\n\nRegarding claim 2)\nThis claim is supported by the results in Table 1 & 2, showing improved accuracy as measured by AUROC and calibration. However, it is unclear to me how strong the baselines are. \n\nRegarding claim 3)\nThe method is claimed to reduce overfitting, and the method does indeed perform well compared to the baseline, but it is unclear whether this is caused by reduced overfitting or something else. Using \"great performance on train set, poor perfomance on val/test set\" as a loosse definition for overfitting, I don't see much evidence for this claim. \n\nRegarding claim 4)\nIt is claimed that the method has a reduced number of tunable hyperparameters. However, the proposed prior has several hyperparameters, e.g. K, c, d_k, e_k. The paper claims that the results are not sensitive to K, but there is no evidence to support this. Moreover, it is also unclear how sensitive the results are to the remaining hyperparameters.", "audience_yn": "Yes", "audience_explanation": "Causality is an important topic and I am sure the a subset of the TMLR audience would be interested in work like this.", "requested_changes": "Overall, the work seems solid, but there are several things in the presentation and the results that makes it harder to appreciate. I'll elaborate on some of the things that would improve the paper in my opinion below.\n\nNotation, symbols, abbreviations are in many places not defined explicitly, e.g.\n- In eq. (1), does A^\\tau refer to the matrix A raised to the power of \\tau or does \\tau index a collection of matrices?\n- Property 3.3 mentions a \"conditional posterior of m\", but it is unclear what it is supposed to be conditioned on. \n- Which parametrization of the gamma distribution is used in eq. (3)?\n- The metric SHD appears only to be defined in the appendix.\n\nThe presentation of the method in sec. 3 appears unclear and unmotivated in several places. E.g.\n-  \"we first present our BVAR model in a general form, which separates the coefficients into binary GC graphs and weight matrices\".\n- \"we propose a new link function named generalized Bernoulli Poisson link (GBPL) that thresholds a random Poisson variable m at V ∈ {1, 2, . . . } to obtain a binary variable b, as one of the key building blocks in our approach.\"\n- Why is it useful to separate graph connections from strength? What is the authors' definition of the a link function? Do they refer to a link function from generalized linear models or something else? What spaces are the link function defined? What problem does it solve?  I can guess some of the answer, but the paper would be improved if the authors would explain, motivate, and justify the construction.\n\nSection 3.4 highlights that the model has fewer hyperparameters, and that a user in essence only has to tune V. The prior has several hyperparameters, e.g. K, c, d_k, e_k and V. How sensitive is the resulting method to use? How does the performance of the method behave as a function of K?\n\nThe proposed prior resembles an Indian Buffet process (IBP) in any ways. It would be great if the authors could compare their proposed prior to the IBP.\n\nThe abstract, introduction, and conclusion highlight interpretability, but it is largely ignored in the rest of the paper. The author should provide some insight into this aspect.\n\n**Experiments**\n\n- Section 5 simply states that the authors focus on AUROC for evaluation, but one needs to go the appendix to figure out what the AUROC is computed for (e.g. link in the graph, where links are estimated using the posterior mean of the Bernoulli variables). Such information needs to be in the main paper in my opinion. \n- Also it is unclear to me what the MSE is measuring. Please elaborate on this.\n- While the authors do show MSE as a function of iterations, it would be useful to know whether the chains have mixed (e.g. via the potential scale reduction factor statistic or similar) and  what the effective samples sizes are etc.\n- It is unclear what the \"uncertainties\" in Table 1&2 represents. It is standard error, standard deviation across the 5 runs or something else?\n- How is the ground truth graph defined for the FMRI dataset?", "broader_impact_concerns": "None", "additional_comments": "The intro says: \"The difficulty of model selection without access to sufficient validation data or ground-truth causal graphs, making hyperparameter tuning unreliable and potentially biased and, ultimately, compromising robustness and interpretability\". Since this is stated as a fact, please support it by citations or own evidence/arguments.\n\nThe intro says: \"Introducing a hierarchical factorized prior over binary Granger causal graphs, which encodes structured sparsity and incorporates inductive bias in the absence of abundant data\". All models/methods have inductive biases, the question is whether they are useful or not. It would be great with a more detail description of the induced inductive bias.\n\nThe Bayesian baseline BVAR is from 1986. Is this really the state of the art for Bayesian approach to Granger causality?", "extra_scores": {}, "extra_text": {}}, {"review_id": "KMqCvcwseo", "reviewer": "Reviewer_X8WU", "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": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper introduces HiBaNG, a hierarchical Bayesian VAR for Granger‑causal discovery in low‑data settings. Coefficients are split into a binary lagged graph ($G_\\tau$) and real‑valued weights ($A_\\tau$). A new Generalized Bernoulli–Poisson Link (GBPL) thresholds a Poisson draw at level (V) to control graph sparsity; this enables a gamma–Poisson factorization prior over $G_\\tau$. Inference uses a conjugate Gibbs sampler with per‑iteration cost $O(N^2)$ and total $O(N^2(V{+}K)\\tau_{\\max}+T\\tau_{\\max}^2)$. On several simulation and real-world tests using climate data, HiBaNG improves AUROC/AUPRC in most cases and shows lower calibration error on most of the settings; a climate (JRA‑55) case study illustrates interpretability.\n\nStrengths\n- Design is reasonable and reproducible: principled uncertainty over edges, simple and transparent inference, good performance in scarce‑data regimes\n- Clear write‑up with hyperparameter grids and clear intuition derivation for the proposed algorithm.\n\nWeaknesses\n- Model selection uses test‑set MSE\n- sensitivity to key hyperparameters is thin\n- noise covariance is diagonal (not realistic)\n- real‑data validation beyond the climate case is limited (i.e. lacks broader real-world impact validations)", "claims_and_evidence_yn": "Yes", "claims_explanation": "- Methods and assumptions are explicit with sufficient references to past works.\n- The Gibbs updates and complexity are derived clearly. \n- Across benchmarks, HiBaNG is typically best or second‑best in AUROC/AUPRC; on FMRI it ranks second (PCMCI+ is first), thus indicating strong support for the statements made in the paper.\n- The experiment setups are well elaborated (e.g. graphical illustrations and tabular results in the appendix, results are computed in multiple runs with mean/variance calculated).", "audience_yn": "Yes", "audience_explanation": "The paper tackles a common pain point \"causal discovery with short time series\" using a compact Bayesian model that returns calibrated edge probabilities. The GBPL‑based prior and the separation of $G_\\tau$ from $A_\\tau$ are likely to interest both causal‑discovery and time‑series researchers.", "requested_changes": "- If possible, please replace test‑set tuning with a validation split or rolling origin CV for all methods\n- Report R‑hat/ESS or multi‑chain traces for key variables and include reliability diagrams for edge probabilities (beyond CE)\n- Vary $V$, $K$, and $\\tau_{\\max}$ for additional ablation studies. You can also discuss the effect of the diagonal covariance assumption (e.g. compare with full covariance)\n- Fix an inconsistency: the text/Appendix tune $V\\in \\{1,2,3\\}$, but Table 3 shows $V=\\{1,3,5\\}$.", "broader_impact_concerns": "The paper has a Broader Impact section. \nTwo items need clearer emphasis: \n- Granger != interventional causality; results hinge on stationarity and unconfounded lags, so edges should be vetted by domain experts, especially in high‑stakes domains. \n- For human or sensitive data, note bias/privacy risks and basic mitigations.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)