DoFlow: Flow-based Generative Models for Interventional and Counterfactual Forecasting on Time Series
Authors: Dongze Wu
Feng Qiu
Yao Xie
Venue: ICLR 2026 Poster (2026)
Keywords: Time Series
Causal Inference
Generative Models
Flow Matching
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)

Abstract: Time-series forecasting increasingly demands not only accurate observational predictions but also causal forecasting under interventional and counterfactual queries in multivariate systems. We present DoFlow, a flow-based generative model defined over a causal Directed Acyclic Graph (DAG) that delivers coherent observational and interventional predictions, as well as counterfactuals through the natural encoding–decoding mechanism of continuous normalizing flows (CNFs). We also provide a supporting counterfactual recovery theory under certain assumptions.  Beyond forecasting, DoFlow provides explicit likelihoods of future trajectories, enabling principled anomaly detection. Experiments on synthetic datasets with various causal DAG structures and real-world hydropower and cancer-treatment time series show that DoFlow achieves accurate system-wide observational forecasting, enables causal forecasting over interventional and counterfactual queries, and effectively detects anomalies. This work contributes to the broader goal of unifying causal reasoning and generative modeling for complex dynamical systems.

Published as a conference paper at ICLR 2026
DOFLOW: FLOW-BASED GENERATIVE MODELS FOR IN- TERVENTIONAL AND COUNTERFACTUAL FORECASTING ON TIME SERIES
Dongze Wu Yao Xie
H. Milton Stewart School of Industrial and Systems Engineering Georgia Institute of Technology
Atlanta, GA 30332, USA dwu38 88 @ga weweeeechcuu.eeeeeeeeeee yao. xxe.iielise www.eeiisyeggttech..e.
Feng Qiu
Northwestem- Argonne Institute for Scientific and Engineering Excellence Evanston, IL 60208, USA Eglu@@ootthwesteen.dde
ABSTRACT
Time-series forecasting increasingly demands not only accurate observational pre- dictions but also causal forecasting under interventional and counterfactual queries in multivariate systems. We present DoFlow, flow-based generative model defined over a causal Directed Acyelic Graph (DAG) that delivers coherent observational and interventional predictions, as well as counterfactuals through the natural en- coding-decoding mechanism of continuous normalizing flows (CNFs). We also provide supporting counterfactual recovery theory under certain assumptions. Beyond forecasting, DoFlow provides explicit likelihoods of future trajectories, enabling principled anomaly detection. Experiments synthetic datasets with various causal DAG structures and real-world hydropower and time series show that DoFlow achieves accurate system-wie wide observational fore- casting, enables causal forecasting over interventional and counterfactual queries, and effectively detects anomalies. This work contributes the broader goal of unifying causal reasoning and generative modeling for complex dynamical systems.
INTRODUCTION
Forecasting the evolution of multivariate time series central problem in statistics and machine learning, with extensive development across both classical and modern approaches. In the standard setting, a model observes past values and predicts what future values are likely to be under the same conditions that generated the past data. These models are purely observational: they learn correlations from past behavior and extrapolate them into the future.
However, many real-world applications require more than observational forecasting. Practitioners often need to answer "what if" questions about interventions, i.e., how system would evolve certain control variables were set differently. In this work, we focus on two types of causal queries:
An interventional query asks: "How will the forecast change under a planned modification of certain variables?" Here, the goal is to predict the future trajectory under a chosen sequence of actions. For example, hydropower time series, turbine-control signals are upstream causes of downstream signals such as power output (Klein et al. (0221). Operators may specify hypothetical turbine-control plan and ask how all downstream signals will evolve. This contrasts with standard observational forecasters, which, given fixed past context, produce fixed forecast. They do not allow one to modify future control plans and inspect how changes propagate through the system. counterfactual query instead asks: "What would the system trajectories have looked like had we intervened differently?" Here, have already observed factual trajectory and ask

---

Published as a conference paper at ICLR 2026
how the same trajectory would have changed under alternative intervention. For example, in healthcare (Cheng et al.; [2020] Jacob Rodrigues et al. 2020). we may observe a patient's 's treatment and outcome trajectory, and then ask whether this particular patient's outcome trajectory would have been better or worse under a different dosing schedule. Outcomes depend not only on dosing but also on patient-specific factors that are not directly observed (e.g., baseline health), yet these factors are implicitly reflected in the factual trajectory. By conditioning on the factual trajectory, we can infer such unobserved factors and then re-simulate how this same patient's 's trajectory would have evolved under the alternative dosing plan.
However, most modern forecasters are observational: they excel at correlational prediction but cannot answer causal queries. Also, to our knowledge, there exists no general framework for counterfactual time-series forecasting yet. Bridging these regimes calls for model that both causally structured and generative, capable of generating consistent and system-wide trajectories under causal queries.
Contributions. In this paper, we address this gap by proposing DoFlow, a generative model based on Continuous Normalizing Flows (CNFs) that explicitly embeds causal directed acyclic graph (DAG) structure. Leveraging the invertibility of CNFs and the temporal conditioning in Neural ODEs, DoFlow provides a unified framework for observational, interventional, and counterfactual forecasting. Beyond forecasting, it also yields explicit likelihoods for future trajectories, enabling principled anomaly detection. By unifying causal reasoning with deep generative modeling, DoFlow advances toward trustworthy inference and decision support in complex dynamical systems.
1.1 RELATED WORK
To situate our work, we briefly review related lines of research spanning time-series forecasting and causal generative modeling.
Time-series forecasting. Modern approaches can be broadly categorized into four families; ii)
classical statistical models such ARIMA, state-space models, and vector autoregression (Young)
and Shellswell 1972, Harvey] 1990; Hyndman et al.] 2008 Zivot and Wang 2006 Scott and
Varian [2144] Taylor and Letham [2018); (ii) deep sequence models including RNNs/LSTT and
attention
attention-augmeened augmented variants (Medsker et al. 2001: Graves 2012 Qin et al.) 2017; Lndemaan
al.
2021);
I) (iii) transformer-bsed architectures offering scalability on long sequences (Vaswani
et al.
2017 Zhou et al., 2021 Wu et al. 2021 Nie et al. 2022); and (iv) recent flow and diffusion
models that directly generate trajectories, including autoregressive (Rasul et al. 2020 2021 Gao et al. 2025). graph-augmented (Dai and Chen 2022), and horizon-wide designs (Alcaraz and Strodthoff 2022, Shen and K Kwok] 2023). These models typically condition on learned context representations and use expressive generative backbones to produce future forecasts, but they remain observational, lacking the ability to simulate interventional and counterfactual trajectories.
Causal generative modeling. For static data, interventional and counterfactual queries are tradi- tionally addressed via structural equation models and do-calculus. Recent generative extensions include graph causal encoders (Schölkopf et al., 2021 (anchezz-arttn et al.
2022] Rahman and Kocaoglu 2024), diffusion and flow-based models (Khemakhem al. 202 Javaloy et al. [023] Chao et al. [2224] Wu et al. 2025), and model-agnostic counterfactual generators (Karimi et al. 2020 Dasgupta et al. 2025). ). However, these methods largely target non-temee non-temporal settings and do not capture causal dependencies over time. A related work studies counterfactual dynamic forecasting (Liu et al. 2023) in constrained physical systems, with interventions only on the initial constrained state. By contrast, DoFlow extends causal generative modeling to general time-indexed DAGs, enabling coherent interventional and counterfactual forecasting across dynamic systems.
We briefly situate DoFlow relative to two complementary lines a practitioner may combine with our method: (i Causal effect estimation-quaaat quantifying how external actions change short-term expected outcomes, and (ii) Causal discovery a causal DAG from observational data.
Causal effects on time series. Prior research on time-series interventions has primarily focused on estimating treatment effects under external actions, using variational Bayesian models, and
classical time-series approaches like ARIMA (Brodersen et al. 2015; Scott and Varian 2014
Linden [0155] Bernal et al.] 2017: Bica et al., 2020b al Menchetti et al.) |021 Wijn al. 2022
Grecov et al.) 2022; Giudice et al. 2022; Runge et al., 2223
and Rostami- Tabar 2025).
These studies typically model how a discrete action variable A, E (0...., J-1) for example,
2

---

Published as a conference paper at ICLR 2026
representing treatment versus control different policy levels affects future outcomes of the series. These are often formalized as тę = | A = 4 = k], where Y, denotes the outcome and A, the time-varying action. Methodologically, modern approaches estimate the
expected outcome under alternative actions, often through conditional generative or variational
representation-leamiie models, to compute average individualized treatment effects (Lim) 2018
Liu et al. 2020 Li et al. 2021 Cao et al. 2023; Wu et al. 2024b Cinquini et al.) 2025).
Causal discovery. Inferring causal directed acyelic graphs (DAGs) from observational time series has also been actively explored. Existing approaches include statistical dependence- tests (Melkas al. 2021 Bussmann et al. 2021], owtimiiationn
formulations (Pamfil ct al., 2020), and deep learning- based frameworks (Tank ct al. 2021; Yin and Baruccal 2022] Zhong et al.,
2023 Cheng al. Faruque et al. 2024 Cheng al.) 2024), which collectively advance the
recovery of temporal causal structures under various assumptions.
However, existing causal effects on time-series methods foeus mainly on discrete, fixed-time actions and estimate short-term expected outcome differences. Moreover, to our knowledge, existing work has modeled the counterfactual trajectory, despite its importance for reliable decision making. Therefore, there remains need for generative model that explicitly embeds causal DAG in time-series, supports interventions on individual continuous variables at arbitrary times, and yields coherent interventional and counterfactual predictions across the system.
2
TIME-CONDITIONED FLOW ON DAG
2.1 SETTINGS AND GOALS
We consider K-dimensional multivariate time series evolving over causal directed acyclic graph (DAG) with nodes {l...., ...,,, K} in a topologically sorted order. We denote by X;, E R the value of node i time-series step t. Let X == (X.,,...., x denote the collection of all nodes time t, and let X pa = {X,,* E } denote the values of node i's parents at time t.
The parent set pa(i) encodes the structural dependencies that remain fixed-for instance, in physical system, upstream components serve as permanent causal parents of downstream ones. At each step t, the value X depends on its past trajectory Xi and the past trajectories of its parents Xpa(1),-;
, (X...* 8 < f), Xpa(i),-- {x,.,, j E pa(i), < t}-
(1)
We assume no tiptst- causal effects, meaning that all causal influences occur with at least one time-step lag. We formalize the dynamics by the structural causal model (SCM):
X1, == , Xpa(s).-, U..))
where Uz, is an exogenous noise independent across nodes and time. Here, fa is a causal mechanism that takes as inputs the histories of and its parents, plus an exogenous noise.
For training and evaluation, each sequence (X....., X is divided into a context window {X....., } used for conditioning and a forecasting window {Xr+1, Xr} used for prediction. The context window consists of observed past values and is never intervened on. Interventions, if present, are applied within the forecasting window and influence all subsequent steps that we aim to predict. We use uppercase for random variables X, and lowercase for realizations X,.
Goals. Our goal to develop framework for causal time-series prediction that answers:
Interventional forecasting. We denote the intervention schedule by I S [K] x {r 1,..., T}. where each () E I specifies an intervention node at time 1. For each (i, L) E EZ, the variable Each intervention
Xz, intervened to the value Ti, forming the collection 7 == replaces the system's natural evolution with fixed value do(X,, i,t == N,,), collectively written as do Xx == 7.). This defines the interventional distribution:
+X,+ X1,7, do(Xz 7z)).
When I = 0, the model reduces to standard observational a forecasting.
Counterfactual forecasting. Given factual trajectory x we aim to answer: "What would this specific sequence have looked like if variables Xz had instead been set to zz?" This induces the counterfactual distribution:
p(X Xi:r, x *T1TT do(Xz = =xx))-
3

---

Published as a conference paper at ICLR 2026
2.2 CONTINUOUS NORMALIZING FLOW FOUNDATIONS
We first introduce the continuous normalizing flow (CNF)) foundation of DoFlow. It autoregressively generates future values, with each step conditioned on past temporal information. For each node in the causal DAG, we learn a lime-conditioned CNF that models the conditional distribution of node i each time step t, given its own and its parents' past temporal context. This context is summarized by hidden state H1, which aggregates information from the historical sequences xx and X pa(t), In practice, H obtained using RNN-based encoders, as detailed in Section 2.3
At each time-series step t, a CNF defines a continuous transformation between the target distribution (at s=0, corresponding to the target data X,,) and a base distribution (at s=1, typically (N 1)) through neural ordinary differential equation (ODE) (Chen et al. 2018). Here, denotes the continuous ODE time, which is distinct from the time-series index
d = (x( 8; H1,-1), [0
(5)
ds
where the velocity field e is parameterized by neural network. To predict T₁,, we sample z (1,r(1)~~W)) 1) and integrate the ODE from & = 1 to s = 0. We train the velocity field Us using Conditional Flow Matching (CFM) loss (Lipman et al. 2023). which directly regresses o onto an analytically defined reference velocity field. Specifically, reference path é interpolates between a data sample di, and a base sample = (typically using linear interpolation), whose derivative 0, defines the reference velocity. The model minimizes an Le loss between o and 0, allowing efficient training while preserving the invertible mapping from base distribution to data distribution.
2.3 TIME-CONDITIONED CONTINUOUS NORMALIZING FLOW
RNN to summarize past histories. We learn continuous normalizing flow (CNF) for each node i to autoregressively predict Xi, for / E {T + 1, . T}, conditioning each step on the node's past history and its parents' past histories. An RNN (can be either LSTM or GRU (Dey and Salem 2017 is employed to summarize histories of node via a recurrent state ha.r
h RNN( hi.
Therefore, for each node at time t, DoFlow is conditioned on the concatenated hidden states:
Hi == concat (hi,c-1; hpp()))-2)),
where hpp(0))-1 == (haa-1 1 : E pa))))
At forecasting time, we use H.,, (for observational or interventional forecasting) or HC (for counterfactual forecasting) to denote the hidden states updated from the model predicted values. Importantly, since f1,,-1 (or H,,,;,) is autoregressively updated from the model's generated or intervened values, it carries information about any past interventions applied node or its parents. Besides, we denote H.,, as the factual hidden state computed from the observed factual trajectory.
Time-conditioned CNE
We train a separate CNF for each node that are shared across time-series steps f, conditioning each step on H.,--1- The Neural ODE of the conditional CNF is defined as:
do),e()) = (x),,(()), 8; H;,-1)),
(8)
ds
s E (1+1,,,, T). which connects samples from the time-series distribution at time (at = 0) with the base distribution (0, 1) (at = 1).
Training loss. At the start of the forecasting window, H is initialized from the context sequence (X,,1., I:r) For later steps t 2 + 1, the hidden states H+,, are updated autoregressively during training using the observed values (Xa,, X Xpa()).) The training loss is calculated over the entire forecasting window. For each t E {r + I,..., TT), we define the reference path 6 as straight- interpolation between the training sample zi,e and Gaussian base sample 2 (0, I):
ol (0 2; = (1 (-s) 5) Z4,e + s2, 0, 2; = Fi,t-
(9)
with U100 1).
4

---

Published as a conference paper at ICLR 2026
Here, 0,6 represents the reference velocity field, which the model's learned velocity field i, aims to approximate. Training thus minimizes the squared Lg distance between v and 0.Ø. Therefore, the training loss of the flow using conditional flow matching (Lipman et al. 2023) becomes:
K
T
(CF = Exirrpx
Σ
K(T
(10)
E
(0(11,));))) s; Hi,t-1) 0, 10% ==,2 (0000)
where I encompasses both the parameters of the velocity field 82 and the RNN parameter.
After training, we obtain a velocity field v; for each node on the causal DAG. We next define the forward and reverse processes of the Neural ODE.
Forward process. We treat the forward process as an encoding operation, denoted by the function Z: (X,;; H1-1). Given the velocity field 2; the forward process pushes an observed factual
-F outcome F (at s = 0) to latent embedding Zi, (at s = 1), conditioned on factual hidden state H. 1,t-1 that summarizes the past observed factual data. Formally,
Do( H.,,- 2T,,, + //'00, (e((,,e()) s; H ds,
(11)
with (0, x.r Through this process, 2,, encodes the information of the factual observation s.,,-
Reverse process. We treat the reverse process as a decoding operation. is initialized differently for interventional (Section and counterfactual prediction (Section (32).. In general, to predict the value of node i at time t, given a latent representation 2i,t, which can be either sampled from A(00 1) or obtained by encoding factual sample, the reverse process is defined as:
i (z,, Ha,-1) 2,,1 -f, ua (Z(,))()), s; H1.t-1)d ds,
(12)
with x;,( (1) zi,x- Here, He. is autoregressively updated using the predicted values ăे. and pa(()) and serves as the conditioning state for the next time step.
INTERVENTIONAL AND COUNTERFACTUAL PREDICTIONS
3.1 OBSERVATIONAL AND INTERVENTIONAL PREDICTION
At inference time, we forecast each node's value one step ahead using the reverse process condi- tioned on the latest hidden states. Assume that our intervention schedule is I with intervened values {Tr.). J. and the purely observational case
Algorithm 1: Time Series Observational/In- terventional Forecasting
1: Input: Context window { ME intervention schedule I with values (m.))
2: Initialize hidden states fl.r with
is given by I = 0.
3: for t 7 (pa,):r;=pa(i)1:> for all i 1,..,,
do for = 1,..., K do (topological order)
The forecasting at time step is proceeded in topologically sorted order, meaning that parent nodes forecasted first, followed by their chil- dren. For node intervened at time with value Na,t, the forecast is fixed to the intervened value,
4:
S:
if (i, I then
6:
54,t +- 7/.t
7:
else
8:
Sample 2 NN(0.)) 1)
i.e. L.. + Ti For non-intervened nodes, the forecast f.,, is generated by the flow model us- ing reverse process (Eq, [2). conditioned on its hidden state H,,-,, which is updated using the previous forecasts -, and & The overall procedure for observational and interven- tional forecasting over the time series is summa- rized in Algorithm [] An illustrative figure is shown Panel B of Figure]]
9:
f.tt $ (zi,e; Alat-1) (Eq.
10:
end if
h4, + (NN(
H. (hr,er
13: end for
end for 15: Output: ( 1=1+1....
5

---

Published as a conference paper at ICLR 2026
RNN State
Observational Interventional
C. Counterfactual
( (-:8ppp) M
Ou Me
Nyye RNN
Zu
F 1
Nit-1
Figure I: (A) RNN State Update. (B) Observationa//nteeeennioon Forecasting. Forccasts are generated by decoding from latent 2,, ~ N (0, 1). conditioned on H updated with the last predicted (áx,e-1, pps(t)))-1) (C) Counterfactual Forecasting. factual observation x encoded with factual state H H into = then decoded under the counterfactual state H to yield a. Factual states H are updated from observed (pa(())--)) while counterfactual states f are updated from the previously predicted (di,e-1 CPP pa(t).t-1)-
factual (2 (2, oo
3.2 TIME SERIES FORECASTING FOR COUNTERFACTUAL QUERIES
Counterfactual forecasting follows the standard Abduction-actions preddctionn procedure. Given factual trajectory {-f, + /==1 we proceed as follows: (i) abduction we infer latent variables by encoding each observed factual value into its latent representation = through the forward process conditioned on a factual hidden state H that summarizes past factual observations; (ii) action apply the specified intervention schedule I; and (iii) prediction
[1,
generate the counterfactual trajectory CCP through the reverse process 12). starting from the abducted latent
representations =
Specifically, we first compute the factual hidden states [ for each from the observed 1
context f 1KK and the observed factual trajectory (-f.-1-T)) These factual states are used only for encoding factual values into their latent representations, which serve the starting point for counterfactual decoding. At each time-series step I, nodes are predicted in topological order. For ti I, the counterfactual forecast is CC == Y.s. and its RNN state hs is updated accordingly. For non-intervened nodes, the factual value =F,, is first encoded as 2 = (xF,;; (f(,,-). then decoded under the counterfactual hidden state o 2, (zT,;; 11, followed by updating CCF The overall counterfactual generation procedure is summarized in Algorithm and illustrated in Panel of Figure
Algorithm 2: Counterfactual Time Series Generation
1:
Input: Context window { 1.K factual sample {2f 1:r 1.K intervention schedule I with values (Tre)
2: Obtain factual hidden states (H, |T, for each i from context { K and observed factual (c],r41.7)4
3: Initialize counterfactual hidden states H with context (24,177 (pa()) all 1, 1....,K
4: for f= + to do
5: for i 1, ...... K do (topological order)
6: if (i, E I then
7:
i +- Ys,s
8:
else
9:
=tt + o H.,..)
(Eq.
11
Abduction)
10:
CCF - ( 1iFE
[Eq. 12 Action-Prediction)
end
12:
har + RNN 1112-1)
13:
HCF (hv,t, ()pa(().))
14:
end for
15:
end for
16:
Output: (e t.r

---

Published as a conference paper at ICLR 2026
3.3 ADDITIONAL PROPERTY: LIKELIHOOD BASED ANOMALY DETECTION
Another advantage of our framework is that it assigns an explicit log-density to its predicted future trajectory over the forecasting window. For single node (index i omitted for clarity), we denote this conditional density by p8, 08.X--1r (- Hr), where H, the context state. In DoFlow, 4:. are base samples used to generate the forecast 17+1TT. The explicit form of this learned density, mapping from latent 21+ predicted 17+1- is given as:
Proposition 3.1. Given base samples 2T-1 q()), the log- log-density of the generated time series obtained via the continuous normalizing flow is:
T
[0g0(2)) - S- vo(2)(6), H.-1) ds
log 11: (1T1IT) H.)
Σ
(13)
where 2, = (2 H1,1), t {T T}.
Anomaly Detection. Since anomalies deviate substantially from normal patterns, expect the model to assign lower density to its generated forecast trajectories when conditioned on anomalous contexts.
4 THEORETICAL PROPERTIES
Our analysis assumes that the underlying SCM on the given DAG satisfies the conditions in Assump- tion [on[B.](ssee Appendix B. which are standard and fundamental in causal inference (Pearl 2009)
Hernán and Robins 2010; Imbens and Rubin 2015 Javaloy et al. 2023 Chao et al. 2024). These assumptions ensure that the interventional and counterfactual distributions we study are well-defined.
In this section, we present theoretical result on the counterfactual recovery of our algorithm under certain assumptions. We fix a node and omit the index for simplicity. Recall that the encoding function is defined as Dg X x H Z and the decoding function 2 Z H - X. corresponding to (11 and (12) respectively. Also, the underlying structural causal model is given by X, := f(X(_ --, pa,t- Ut). In the following, present supporting result on the counterfactual recovery properties DoFlow. We begin by introducing the following assumptions.
Assumption 4.1.
(A1) U2 IL "Ki_,Xpatt-
(A2) The structural causal equation [(, is strictly monotone and continuous in U.-
(A3) For the encoded latent variable Ze = (xx;; H,-1), the conditional distribution satisfies (e( H1-1) (Z)).
Remark 4.2. Each node X, E IR, so (A2) is assumed in the univariate case. It is automatically U,. Under certain identifiability
satisfied under additive SCMs, i.c. X, = (X conditions, it also holds for non nin-linear models (Zhang and Hyvarinen 2012 Strobl and Lasko 2023). Under (A2), for any fixed support of U,. One may notice that (A1)-(A2) mirror those in Bijective Generation Mechanisms (BGM) (Nasr- Wasr-EEfahany et al. 2023). which establish model-agnostic identifiability with an additional assumption on distribution matching. In contrast, our Corollary 4.5 provides recovery for our CNF under (A1)-(33) without requiring distribution matching; see Appendix B.2
the map u f(X, X pm.. is a bijection on the
4.5 model-specific, pointwise
Note that (A3) implies that the encoded Z: &b (X:: Ht-1) is statistically independent of H and equivalently of (X,-,Xpa Xpa.c-), ), in distribution. In the infinite-date limit with exact training, the continuous normalizing flow maps every Xt, conditioned on any fixed (X,_.Xpa,-- ), to the same base distribution ((Z)) = (0, 1), so that (A3) holds exactly.
Under Assumption 4.1
present the first result in this paper:
Proposition 4.3 (Encoded as a function of the exogenous noise U.). Let Assumption 4.1 hold. Without loss of generality, suppose the exogenous noise U, Unif(o, 1). At each time I, the observed variable is generated by the structural causal model X, = f(..,, Xpa,t_ U.), and that the flow encoder produces Z, 46 (Xx;; He-1). Then there exists continuously differentiable bijection g :u Z. functionally invariant to H such that,
Zi = (X Ho-1) = 2 (f Xpa.t, U,); H,-1) = g(Ui) a.s.
(14)
7

---

Published as a conference paper at ICLR 2026
Remark 4.4. Proposition 4.3 states that Z, is function of the exogenous noise Ut, and for con- venience assumes U, ~ Unif [0. 1]. This distributional assumption can be relaxed to other noise distributions. For example, if Z NO 1) with CDF F. then U F(Z) Unif[0, 1], and any assignment U) can equivalently be written as f(,2) ((,)()))
Following Proposition 4.3 we state a counterfactual recovery result under monotone SCMs, Given an intervention schedule I with values {ye), the true counterfactual process is defined recursively as:
E
X CF
MXX pa.t-, Ui), otherwise,
(15)
where U, is abducted from the factual sample XF. Also, X and X pa./-1 denote the model predicted counterfactual values by Algorithm 2 We can now state the following result:
Corollary 4.5 (Counterfactuaal recovery). Let Assumption 4.1 hold. Consider a factual sample generated by the structural causal model X = -X Uz), and let its encoded latent be z = Lo (Xf; Hf-1). At time step t, we obtain the counterfactual hidden state H from the predicted counterfactual past (8 _x pa,t-) The decoder recovers the true counterfactual at time step t almost surely:
x = =6/ ( H = =X
(16)
EXPERIMENTS
In this section, we evaluate DoFlow on observational forecasting, interventional and counterfactual forecasting, and anomaly detection, on both synthetic and real-world datasets. Our code is publicly com http://www...coetttttssss/ooo//a
available at
Sactooloooo__sssseseee Time Trime_Series Series
In real-world, the ground-trum counterfactual is never observable, and the ground-trun interventional only observable if interventions are actively conducted under correct causal DAG. Therefore, we rely on synthetic experiments to obtain quantitative performance metrics. For real-world evaluation, we assess our model on interventional queries and anomaly detection using hydropower datasets from Argonne National Laboratory. and on interventional treatment effect estimation using the cancer-treatment dataset from Bica et al. (2020a).
5.1 SYNTHETIC DATA EXPERIMENTS
We consider both the additive noise model, i.c., x U(,))) =f ff - x UKe, which satisfies Assumption 4.1 and supports the counterfactual recovery result, as well as more general non-linear and non-additive (NLNA) cases of f: ( p robustness beyond the scope of our counterfactual recovery result.
U,,,) to test the model's
We evaluate the model across multiple structurally diverse causal DAGs using Root Mean Squared Error (RMSE), Maximum Mean Discrepancy (MMD), and Continuous Ranked Probability Score (CRPS). RMSE is reported for all settings. MMD and CRPS are computed only for observational and interventional forecasting, since counterfactuals in DoFlow yield a single deterministic trajectory. Detailed simulation setups and metric definitions are provided in Appendix [
For baseline comparisons in observational forecasting, we consider: a pure RNN-based method, the Gated Recurrent Unit (GRU) (Che al. 2018); ; transformer-based methods, including the Temporal Fusion Transformer (TFT) (Lim et al.
2021 and the Time Series Dense Encoder (TiDE) (Das et al., 2024); and an all-MLP method with specialized contextual mixing structure, the Time Series Mixer (TSMixer) (F.kambaamm et al. 2023). We also compare against probabilistic methods, including DeepVAR (Salinas et al. 2019) a deep RNN-based model with multivariate Gaussian outputs, and convex deep neural network that learns multivariate quantile functions. We extend our acknowled Lekkowleegemmntss to the Python packages Darts (Herzen et al. 2022) and Gluon (Alexandrov al. 2020). which we use to directly test several modern baselines in this paper.
MQF2 (Kan et al.) 2022)
To our knowledge, comparatively few works tackle interventional or counterfactual time-series forecasting causal DAG. Nevertheless, can adapt strong observational forecasters for interven- tional task by training a separate model for each node in the causal DAG, conditioned on its parents; Appendix E. However, counterfactual generation remains challenging for these baselines.
8

---

Published as a conference paper at ICLR 2026
Table 1: RMSE for observational, interventional, and counterfactual time series forecasting across causal structures: Tree, Diamond, and FC-Layer. Results for the Chain structure are provided in Table 6 (Appendix). MMD and CRPS results are provided in Table and Table (Appendix), respectively. Reported values are averaged over 50 test batches, each containing 128 test series.
Thee
FoLaye
Adiliin
NLSA
Addine
NENA
Addito
NENA
Isl.
Cb..
Is.
Ot.
Obh
Ct..
CF
(154..
0.33;
0.77
8.31
0.58. 0.13, 1039
0:31.
0.10.0
0.36.
cau
1001
10
0372.s
0722.11
1:35.
NA.
1.
1
109,
1000g
1:22,
L:3
108;
0835
10.0
1801.
NA.
85.
NA
AR22
122
13
035
1.7,
1.5aaz
NA.
1887au
N4.
Tables (main) and 5 (Appendix) report the RMSE, MMD, and CRPS results for observational, interventional, and counterfactual time-series forecasting across multiple causal structures (see Ap- pendix D). DoFlow consistently delivers strong performance standard observational forecasting and interventional forecasting compared with our adapted baselines, and uniquely supports counter- factual forecasting. We present visual results for both interventional and counterfactual forecasting in Figure 2 (main text), and in Figures 8 and (Appendix).
Figure 2: Left: "Layer" interventional forecasting results. Nodes X1,t, X2.s, and X3.x are intervened. DoFlow provides 50% and 90% prediction intervals; orange lines indicate the true interventional future. Right: "Tree" counterfactual forecasting results. Node X1, is intervened. DoFlow provides a single forecast in green; the orange lines indicate the true counterfactual future.
5.2 REAL APPLICATION: HYDROPOWER SYSTEM
We evaluate DoFlow on real-world hydropower time-series data from Argonne National Laboratory. In this system, water drives turbine that powers generator and passes through transformer before reaching the grid, with control systems monitoring the process. Signals such as water flow, vibration. and electric current are recorded from each component, forming natural DAG shown in Fig. 10
Our evaluation focuses on two tasks: (1) whether DoFlow can accurately forecast each component's time series under interventional queries; and (2) whether DoFlow can accurately detect power outages in advance using log-density, as discussed in Section
Figure 3 illustrates interventional forecasting in the hydropower system during true power out- age. At this point, the turbine signals (X1,,X break down, leading the entire system into forced outage. Notably, DoFlow successfully predicts the characteristic "spikes" in the generator sig- nals (Xx,, and Xs,)) that follow the turbine fail- ure, demonstrating that our DAG-based approach captures the complex causal relationships across system. Quantitative metrics, together with the anomaly detection visualizations from the second task, are provided the Appendix F3
Hydropower System
Interventional Forecasting
(ntervenel)
Innarschedl
Figure 3: Hydropower Interventional.

---

Published as a conference paper at ICLR 2026
5.3 REAL APPLICATION: CANCER TREATMENT OUTCOMES
Our DoFlow can naturally model interventional forecasting for causal treatment effect estimation in time series. The causal effects on time series were thoroughly introduced in Section [.] We evaluate DoFlow on the auncer-treatment benchmark of Bica et al. (2020a), which contains daily patient-level tumor volumes (outcome) and administered therapies (actions). The action variables X,-1 (x 1==1 (chemotherapy assignment, radiotherapy assignment, chemotherapy dosage, and radiotherapy dosage) act as causal parents of the outcome Y, (tumor volume) on the causal DAG, For each test patient, we use the first 55 days as observational context and perform interventional rollouts for the next 5 days under ten distinct treatment options, yielding predicted outcomes Y0+1:4+5- The normalized RMSE between the true Y+14: and predicted Y4+1.045 is reported in Appendix F.4
The detailed settings and results are provided in Appendix F.4
Notably, DoFlow achieves substantial improvements over prior baselines on this causal treatment effect task, as summarized in Table
5.4 DISCUSSIONS AND LIMITATIONS
Training Efficiency. DoFlow trains separate flow per node in the causal DAG, but each network can be much shallower than transformer models. As shown in Tables 2 and its total model size remains modest, and its training and sampling times are comparable to modern baselines.
Interventions. We adapt modern baselines for conditional intervention generation, but DoFlow performs notably better: its RNN- flow design jointly encodes causal histories and propagates interventions through recurrent states, while the flow backbone enables coherent decoding.
Counterfactuals. As discussed in the Introduction, DoFlow among the early approaches for generating counterfactual time-series and causal discovery on time series. Such counterfactual generation is crucial for modeling component relationships and enabling post-hoc analyses in crucial areas.
causal DAG, complementing work on treatment effects
Limitations Known DAG. DoFlow currently assumes a known and correctly specified causal DAG, as well as causal sufficiency (see Assumption B.1. In practice, the DAG may come from known physical relationships or from upstream causal discovery methods as discussed in Section
I.1 A natural direction for future work is to couple DoFlow with time-series causal discovery or
deconfounding causal generative models (Almodóvar et al. 2025 Xia al. 20222, allowing the model to operate learned temporal causal structures for more reliable causal reasoning.
CONCLUSIONS
We introduced DoFlow, a flow-based generative framework for causal time-series forecasting that unifies observational, interventional, and counterfactual queries on DAG-structured causal systems. Moreover, we provide a supporting counterfactual recovery result for DoFlow (Corollary 4.5) under certain assumptions. Experiments various synthetic DAGs and real applications show strong obser- vational forecasts and effective prediction under causal queries. This framework lays the foundation for domains such energy. healthcare, and other areas where counterfactual and interventional forecasting is critical. The integration of causal structure and generative modeling represents a step toward general theory of inference and control in complex dynamical environments.
ACKNOWLEDGMENT
This work is partially supported by NSF CAREER CCF-1650913, NSF 2015787, CMMI-2112533, DMSS-938106 1M3-1830210, and the Coca-Cola Foundation.
CMMI-
REFERENCES
Juan Miguel Lopez Alcaraz and Nils Strodthoff. Diffusion-basse time series imputation and forecast- ing with structured state space models. arXiv preprint arXiv:2208.09339, 2022.
Alexander Alexandrov, Konstantinos Benidis, Michael Bohlke-Schneider, Valentin Flunkert, Jan Gasthaus, Tim Januschowski, Danielle C. Maddix, Syama Rangapuram, David Salinas, Jasper
10

---

Published as a conference paper at ICLR 2026
Schulz, Lorenzo Stella, Ali Caner Türkmen, and Yuyang Wang. GluonTSS: Probabilistic and Neural Time Series Modeling in Python. Journal of Machine Learning Research, 21(116):1-6 2020. URL
httt:://mmı org/ wttp://ppppppp/////aeaeaae 1tt999888.hhtt. htmm
Alejandro Almodóvar, Adrián Javaloy, Juan Parras, Santiago Zazo, and Isabel Valera. Decaflow: deconfounding causal generative model, arXiv preprint arXiv:2503.155 151 14, 2025.
James Lopez Bernal, Steven Cummins, and Antonio Gasparrini. Interrupted time series regression for the evaluation of public health interventions: a tutorial. International journal of epidemiology, 46(1): (6(())344-3355 2017.
Ioana Bica, Ahmed Alaa, James Jordon, and Mihaela van der Schaar. Estimating counterfactual treatment outcomes over time through adversarially balanced representations. In International Conference on Learning Representations, 2020a.
Ioana Bica, Ahmed Alaa, Charles Lambert, and Mihaela van der Schaar. Time series deconfounder: Estimating treatment effects over time the presence of hidden confounders. International Conference on Machine Learning, 2020b.
Kay H Brodersen, Fabian Gallusser, Jim Koehler, Nicolas Remy, and Steven L Scott. Inferring causal impact using bayesian structural time-series models. 2015.
Bart Bussmann, Jannes Nys, and Steven Latré. Neural additive vector autoregression models for causal discovery time series. International Conference on Discovery Science, pages 446-460. Springer, 2021.
Defu Cao, James Enouen, Yujing Wang, Xiangehen Song, Chuizheng Meng, Hao Niu, and Yan Liu. Estimating treatment effects from irregular time series observations with hidden confounders. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 37, pages 6897-6905, 2023.
Patrick Chao, Patrick Blöbaum, Sapan Kirit Patel, and Shiva Kasiviswanathan. Modeling causal mechanisms with diffusion models for interventional and counterfactual queries. Transactions Machine Learning Research, 2024.
Zhengping Che, Sanjay Purushotham, Kyunghyun Cho, David Sontag, and Yan Liu. Recurrent neural networks for multivariate time series with missing values. Scientific reports, 8(1):6085, 2018.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
Li-Fang Cheng, Bianca Dumitrascu, Michael Zhang, Corey Chivers, Michael Draugelis, Kai Li, and Barbara Engelhardt. Patient-speciffc effects of medication using latent force models with gaussian processes. In International Conference on Artificial Intelligence and Statistics, pages 4045-4055. PMLR, 2020.
Yuxiao Cheng, Ziqian Wang. Tingxiong Xiao, Qin Zhong, Jinli Suo, and Kunlun He. Causaltime: Re- alistically generated time-series for benchmarking of causal discovery. In The Twelfth International Conference on Learning Representations.
Yuxiao Cheng, Lianglog 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, pages 11525-11533, 2024.
Martina Cinquini, Isacco Beretta, Salvatore Ruggieri, and Isabel Valera. A practical approach to causal inference over time. In Proceedings of the AAAI Conference on Artificial Intelligence. volume 39, pages 14832-14839, 2025.
Enyan Dai and Jie Chen. augph-augmented normalizing flows for anomaly detection of multiple time series. arXiv preprint arXiv:2202.07855. 2022.
Abhimanyu Das, Weihao Kong, Andrew Leach, Shaan K Mathur, Rajat Sen, and Rose Yu. Long-tem forecasting with tide: Time-series dense encoder. Transactions on Machine Learning Research, 2024.
11

---

Published as a conference paper at ICLR 2026
Sopam Dasgupta, Sadaf Halim, Joaquín Arias, Elmer Salazar, and Gopal Gupta. Mc3g: Model agnostic causally constrained counterfactual generation. arXiv preprint arXiv:2508. 17221, 2025.
Rahul Dey and Fathi M Salem. Gate-variants of gated recurrent unit (gru) neural networks. In 2017 IEEE 60th international midwest symposium on circuits and systems (MWSCAS), pages 1597-1600.
IEEE, 2017.
Xiaowen Dong. Dorina Thanou, Laura Toni, Michael Bronstein, and Pascal Frossard. Graph signal processing for machine learning: A review and new perspectives. IEEE Signal processing magazine, 37(6):117-127, 2020.
Vijay Ekambaram, Arindam Jati, Nam Nguyen, Phanwadee Sinthong, and Jayant Kalagnanam. Tsmixer: Lightweight mip-mixer model for multivariate time series forecasting. In Proceedings of the 29th ACM SIGKDD conference on knowledge discovery and data mining, pages 459-469, 2023.
Omar Faruque, Sahara Ali, Xue Zheng, and Jianwu Wang. Ts-causalon: Learning temporal causal relations from non-linear non-stationary time series data. arXiv preprint 3rXiv:2404.014666 2024.
Jiaxin Gao, Qinglong Cao, and Yuntian Chen. Auto- moving diffusion models for time series forecasting. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 39, pages 16727-16735, 2025.
Jun Gao, Dongze Wu, Feng Yin, Qinglei Kong, Lexi Xu, and Shuguang Cui. Metaloc: Learning to learn wireless localization, IEEE Journal on Selected Areas in Communications, 4113337881-33877 2023.
Gianluca Giudice, Sara Geneletti, and Konstantinos Kalogeropoulos, Inference on causal effects of interventions time using gaussian processes. arXiv preprint arXiv:2210.028500 2022.
Alex Graves. Long short-term memory. Supervised sequence labelling with recurrent neural networks, pages 37-45, 2012.
Priscila Grecov, Ankitha Nandipura Prasanna, Klaus Ackermann, Sam Campbell, Debbie Scott, Dan Lubman, and Christoph Bergmeir. Probabilistic causal effect estimation with global neural network forecasting models. IEEE Transactions on Neural Networks and Learning Systems, 35(4): 4999-5013, 2022.
Andrew Harvey, Forecasting, structural time series models and the kalman filter. 1990.
Miguel A Heráán and James M Robins. Causal inference, 2010.
Julien Herzen, Francesco Lässig, Samuele Giuliano Piazzetta, Thomas Neuer, Léo Tafti, Guillaume Raille, Tomas Van Pottelbergh, Marek Pasieka, Andrzej Skrodzki, Nicolas Huguenin, et al. Darts: User-friendly modern machine learning for time series. Journal of Machine Learning Research, 23 (124):1-6, 2022.
Rob Hyndman, Anne Koehler, Keith Ord, and Ralph Snyder. Forecasting with exponential smoothing: the state space approach. Springer, 2008
Rob Hyndman and Bahman Rostami-Tabar. Forecasting interrupted time series. Journal of the Operational Research Society, 76(4):790-803, 2025.
Guido W Imbens and Donald B Rubin. Causal inference in statistics, social, and biomedical sciences.
Cambridge university press, 2015.
Mariana Jacob Rodrigues, Octavian Postolache, and Francisco Cercas. Physiological and behavior monitoring systems for smart healthcare environments: A review. Sensors, 20(8):2186, 2020.
Adrián Javaloy, Pablo Sänchez-Martin, and Isabel Valera. Causal normalizing flows: from theory to practice. Advances in Neural Information Processing Systems, 36:58833-58866 2023.
Kelvin Kan, François- François-Xavier Xavier Aubet, Tim Januschowski, Youngsuk Park. Konstantinos Benidis, Lars Ruthotto, and Jan Gasthaus. Multivariate quantile function forecaster. In International Conference Artificial Intelligence and Statistics, pages 10603-10621. PMLR, 2022.
12

---

Published as a conference paper at ICLR 2026
Amir- Karimi, Gilles Barthe, Borja Balle, and Isabel Valera. Model- Model-agnostic counterfactual explanations for consequential decisions. In International conference on artificial intelligence and statistics, pages 895-905. PMLR, 2020.
llyes Khemakhem, Ricardo Monti, Robert Leech, and Aapo Hyvarinen. Causal autoregressive flows. In International conference on artificial intelligence and statistics, pages 3520-3528. PMLR, 2021.
Igor Klein, Stefan Mayr, Ursula Gessner, Andreas Hirner, and Claudia Kuenzer. Water and hy- dropower reservoirs: High temporal resolution time series derived from modis data to characterize seasonality and variability, Remote Sensing of Environment, 13:1112 13207, 2021.
Rui Li. Stephanie Hu, Mingyu Lu, Yuria Utsumi, Prithwish Chakraborty, Daby M Sow, Piyush Madan, Jun Li, Mohamed Ghalwash, Zach Shahn, et al. G-net: recurrent network approach to g-computation for counterfactual prediction under a dynamic treatment regime. In Machine Learning for Health, pages 282-299. PMLR, 2021.
Bryan Lim. Forecasting treatment responses over time using recurrent marginal structural networks. Advances in neural information processing systems, 31, 2018
Bryan Lim, Sercan Ô Ank, Nicolas Loeff, and Tomas Plister. Temporal fusion transformers for interpretable multi-horizon time series forecasting. International journal of forecasting, 37(4): 1748 1764, 2021.
Benjamin Lindemann, Timo Müller, Hannes Vietz, Nasser Jazdi, and Michael Weyrich. A survey on long short-term memory networks for time series prediction. Procedia Cirp, 99,650-655, 2021.
Ariel Linden. Conducting interrupted time-series analysis for single-and multiple-grou, comparisons. The Stata Journal, 152)0080--500 2015.
Yaron Lipman, Ricky TQ Chen, Heli Ben-Hamu, Maximilian Nickel, and Matthew Le. Flow matching for generative modeling. In International Conference on Learning Representations, 2023.
Ruoqi Liu, Changchang Yin, and Ping Zhang. Estimating individual treatment effects with time- varying confounders. In 2020 IEEE international conference data mining (ICDM), pages 382-391. IEEE, 2020.
Yanzhu Liu, Ying Sun, and Joo-Hwee Lim. Counterfactual dynamics forecasting-a new setting of quantitative reasoning. In Proceedings of the AAAI Conference Arlificial Intelligence, volume 37, pages 1764-1771, 2023.
Mohammad Ali Mansournia, Goodarz Danaei, Mohammad Hossein Forouzanfar, Mahmood Mah- moodi, Mohsen Jamali, Nasrin Mansournia, and Kazem Mohammad. Effect of physical activity on functional performance and knee pain in patients with osteoarthritis: analysis with marginal structural models. Epidemiology. 23(4):631-640) 2012.
James E Matheson and Robert L, Winkler. Scoring rules for continuous probability distributions. Management science, 22( 1087-1096, 1976.
Larry Medsker, Lakhmi Jain, al. Recurrent neural networks. Design and applications, 5(64-67):
2, 2001.
Laila Melkas, Rafael Savvides, Suyog H Chandramouli, Jarmo Mäkelä, Tuomo Nieminen, Ivan Mammarella, and Kai Puolamäki. Interactive causal structure discovery in earth system sciences. In The KDD'1 '1 Workshop on Causal Discovery, pages 3-25. PMLR, 2021.
Fiammetta Menchetti, Fabrizio Cipollini, and Fabrizia Mealli. Estimating the causal effect of an intervention in time series setting: the c-arima approach. arXiv preprint arXiv:2103.06744 2021.
Arash Nasr-Esfahany, Mohammad Alizadeh, and Devavrat Shah. Counterfactual identifiability of bijective causal models. In International conference on machine learning, pages 25733-25754 PMLR, 2023.
Venkat Nemani, Luca Biggio, Xun Huan, Zhen Hu, Olga Fink, Anh Tran, Yan Wang, Xiaoge Zhang, and Chao Hu. Uncertainty quantification in machine learning for engineering design and health prognosties: A tutorial. Mechanical Systems and Signal Processing, 205: 10796, 2023.
13

---

Published as a conference paper at ICLR 2026
Yuqi Nie, Nam H Nguyen, Phanwadee Sinthong, and Jayant Kalagnanam. A time series is worth 64 words: Long-term forecasting with transformers. arXiv preprint arXiv:2211.14730, 2022.
Roxana Pamfil, Nisara Sriwattanaworachai, Shaan Desai, Philip Pilgerstorfer. Konstantinos Geor- gatzis, Paul Beaumont, and Bryon Aragam. Dynotears: Structure learning from time-series data. In International Conference on Artificial Intelligence and Statistics, pages 1595-1605. Pmlr, 2020.
Judea Pearl. Causality. Cambridge university press, 2009.
Yao Qin, Dongjin Song, Haifeng Chen, Wei Cheng, Guofei Jiang, and Garrison Cottrell. A dual-stage ceanntion-baeee recurrent neural network for time series prediction. arXiv preprint arXiv:1704.02971, 2017.
Md Musfiqur Rahman and Murat Kocaoglu. Modular learning of deep causal generative models for high-dimensional causal inference. arXiv preprint arXiv:2401.01-23 2024.
Kashif Rasul, Abdul-Saboor Sheikh, Ingmar Schuster, Urs Bergmann, and Roland Vollgraf. Mul- tivariate probabilistic time series forecasting via conditioned normalizing flows. arXiv preprint 20Xx::2002.06103, 2020.
Kashif Rasul, Calvin Seward, Ingmar Schuster, and Roland Vollgraf. Autoregressive denoising diffusion models for multivariate probabilistic time series forecasting. In International conference on machine learning. pages 8857-8868. PMLR, 2021.
Jakob Runge, Andreas Gerhardus, Gherardo Varando, Veronika Eyring, and Gustau Camps-Valls. Valls. Causal inference for time series. Nature Reviews Earth & Environment, 4(7):487-505, 2023.
David Salinas, Michael
Laurent Callot, Roberto Medico, and Jan Gasthaus. High- dimensional multivariate forecasting with low-rank gaussian copula processes. Advances in neural information processing systems, 32, 2019.
Pablo Sánchez-Martin, Miriam Rateike, and Isabel Valera. Vaca: Designing variational graph autoencoders for causal queries. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 36, pages 8159-811 88699,88887188 2022.
Bernhard Schölkopf, Francesco Locatello, Stefan Bauer, Nan Rosemary Ke, Nal Kalchbrenner, Anirudh Goyal, and Yoshua Bengio. Toward causal representation learning. Proceedings of the IEEE, 109(5):612-66)) (2-5))612-3444 634, 2021.
Steven L Scott and Hal R Varian. Predicting the present with bayesian structural time series. International Journal of Mathematical Modelling and Numerical Optimisation, 5(1-2):4-23, 2014.
Lifeng Shen and James Kwok. Non-autoregressive conditional diffusion models for time series prediction. In International Conference Machine Learning, pages 31016-31029. PMLR, 2023.
Eric V Strobl and Thomas A Lasko. Identifying patient-specific root causes with the heteroscedastic noise model. Journal of Computational Science, 72: 72:102099, 2023.
Alex Tank, Ian Covert, Nicholas Foti, Ali Shojaie, and Emily B Fox. Neural granger causality. IEEE Transactions on Pattern Analysis and Machine Intelligence, (00))2267-42999 2021.
Sean J Taylor and Benjamin Letham. Forecasting at scale. The American Statistician, 72(1):37-45. 2018.
Yifeng Tian, Nishant Panda, and Yen Ting Lin. Liouville flow importance sampler. arXiv preprint arXiv:2405.06672, 2024.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones. Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention all you need. Advances in neural information processing systems, 30, 2017.
Stan RW Wijn, Maroeska M Rovers, and Gerjon Hannink. Confounding adjustment methods in longitudinal observational data with a time-varying treatment: a mapping review. BMJ open, 12 (3)e0059977, 2022.
14

---

Published as a conference paper at ICLR 2026
Dongze Wu and Yao Xie. Annealing flow generative models towards sampling high- digh-dimensional and multi-modal distributions. In Forty-second International Conference on Machine Learning, 2025.
Dongze Wu, Jun Gao, and Feng Yin. Bayesian- boosted metaloc: Efficient training and guaranteed generalization for indoor localization. In ICASSP 2024-2024 IEEE International Conference Acoustics, Speech and Signal Processing (ICASSP), pages 6070-6074, IEEE, 2024a.
Dongze Wu, David Inouye, and Yao Xie. Po-flow: Flow-based generative models for sampling potential outcomes and counterfactuals. arXiv preprint arXiv:2505. 16051, 2025.
Haixu Wu, Jiehui Xu, Jianmin Wang, and Mingsheng Long, Autoformer: Decomposition transformers with auto-correlation for long-tern term series forecasting, Advances in neural information processing systems, 34:22419-22458 2021.
Shenghao Wu, Wenbin Zhou, Minshuo Chen, and Shixiang Zhu. Counterfactual generative models for time-varying treatments. In Proceedings of the 30th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, pages 3402-3413, 2024b.
Kevin Xia, Yushu Pan, and Elias Bareinboim. Neural causal models for counterfactual identification and estimation. arXiv preprint arXiv:2210.0003: 2022.
Zexuan Yin and Paolo Barucca. Deep recurrent modelling of granger causality with latent confound- ing. Expert Systems with Applications, 207:118036, 2022.
Young and Stephen Shellswell. Time series analysis, forecasting and control. JEEE Transactions on Automatic Control, 17(2):281-283 1972.
Kun Zhang and Aapo Hyvarinen. On the identifiability of the post-nonlinear causal model. arXiv preprint arXiv: 1205. 2012.
Qin Zhong, Yuxiao Cheng, Zongren Li, Dongjin Wang, Chongyou Rao, Yi Jiang, Lianglong Li, Ziqian Wang, Pan Liu, Yawei Zhao, al. Ultra-efficient causal deep learning for dynamic csa-aki detection using minimal variables. medRxiv, pages 2023-12, 2023.
Haoyi Zhou, Shanghang Zhang, Jieqi Peng. Shuai Zhang, Jianxin Li, Hui Xiong, and Wancai Zhang. Informer: Beyond efficient transformer for long sequence time-series forecasting. Proceedings of the AAAI conference artificial intelligence, volume 35, pages 11 11106-1115. 2021.
Eric Zivot and Jiahui Wang. Vector autoregressive models for multivariate time series. In Modeling financial time series with S-PLUSØ, pages 369-413. Springer, 2006.
A PROOFS
Proposition 3.1 Given base samples 2141:T q(.), the log-density of the generated time series obtained via the continuous normalizing flow is:
[ogg(z)) + S'" V. vo(2r(8); Ar-1) ds
Po,x, 11:0 (Br+IT) M..) = Σ
t=T+1
where G. 2 Hr-1), LE EE +T
Samples from Continuous Normalizing Flows (CNFs) evolve according the following Neural ODE:
Proof.
(s) = v(xr(s), ff,-1), E [0, 1].
(17)
ds
which induces a corresponding evolution of the sample density governed by Liouville's continuity equation (Chen et al. 2018):
Opa(Te, s
lir1))+ (po(ze, Aie ( s; (i-1)) 0.
(18)
Here, po (fx(())) 8) denotes the probability density of the sample (s) ODE-time
15

---

Published as a conference paper at ICLR 2026
Next, we have that the dynamics of the density pa(-) governed by the velocity field va He-1) is given by:
d
Vpa (g
Dspo
log (x)))( H(-1)
(19)
ds
Pe
(xx()).e)
Vpo V (pov)
(by 17 and 18)
(20)
Pø
(2)
Vpa v (Va 1 pa
(21)
Pa
(x)()).8)
(22)
Starting from a base sample 2 ~ q(.) and integrating from & = to = 0, we have:
log (&e (i--1) log ((z)) + L. Vo v(z((5), s; fil-1) ds.
(23)
By summing up all the log-densities within the forecasting window, we obtain:
T
log (2)) + L/' V. us (x((s)), (l--)) ds
log Pe.x ( 11r))
Σ
(24)
1=+1
where ū; = 0g D6 (2:; H1, LE (1,....
Proposition 4.3 [Eccdddd function of the exogenous noise U). Let Assumption 4.1 hold. Without loss of generality, suppose the exogenous noise U, Unif(0, 1]. At each time the observed variable is generated by the structural causal model X, Xpa,t- U.). and that the flow encoder produces Z, = Do (Xe; H2-1)- Then there exists a continuously differentiable bijection u - Z. functionally invariant H such that,
Ze (oo(x:: H2-1)
U.); U,):H,1)) (X(_(XX-,
=g(U)) a.s.
Proof. Fix a node i (index suppressed) and a time t. We write the extended parental state as S. = (X,,, Xpatt) Since the hidden state depends solely on St, we define
qs. (U,) = Z, =
U.): H,-1).
Therefore, our goal becomes proving that Z, qs, (U6) a function invariant S.. By (A3) in Assumption 4.1 we have Ze = qs, (Ut) 1L S₁, and thus:
P2, (2)) pz, (z+).
(25)
Because continuous normalizing flows are invertible, the encoding function Ty X' H - 2 is
invertible with respect to X. In I-D, this implies monotonicity of 4e X, f(X;,; Xpa,t- U.). Without loss of generality,
assume that dy is strictly increasing in X, f. Moreover, by (A2),
since (-, U2) is strictly increasing Ug, follows by the composition rule that qs, (U,) strictly increasing in Ug and hence bijective in [0, 1].
Since Z, = qs. qs. (U+), we may apply change of variables formula:
P (zi) = (u.(,,, de, 9,, (2) =1
dz, 9,, (z),
(26)
where the last equation follows from the uniform distribution of U, and the fact that Pz, 0. The absolute value is dropped because of the (WLOG) assumption that qe, is strictly increasing.
Then because of 25 b. we have that for fixed 2t, d q,, (2, is the same for any pair 8, (p,,, follows that:
90 (xt) = 1
dr - (x)dr = =(xddx)dx
d
c(z)dr Co.
(27)
16

---

Published as a conference paper at ICLR 2026
where c(z))
d dz, As. (2x) is independent of ss, and Cs2 is constant for each St-
By re-inverting (27). we have:
qs. (we) = (90, (goo.)-'(uu) (ue) = (f'o((t) c(x) dr tcs.) [(x))d
""())((()))))) (u))=G GG- Cei)
(28)
:=(())
Since qs, a bijection [0, 1] supp(Z,) Z,) = {2 E R pz, (2)) > 0}, we have:
qs (0) inf supp(Zz), qs. (1) sup supp Zi), for any state st.
(29)
Therefore,
(0) = G (-c,,) = inf supp(Z:).
(30)
The support Z does not depend on S. because of (25). As result, for any state st = (Iz- xpa.t we have that Cs2 == e is constant that does not depend on 8t.
Therefore, we can write
o(U,) G
(U- -c)
9s. (U+), 4ss.
(31)
As result, we conclude that
Ze Do (Xti H1-1) (U,( g(U)) a.s..
(32)
Corollary 4.5 (Counterfactual recovery). Let Assumption 4.1 hold. Consider a factual sample generated by the structural causal model X = f(X,-, X X U2), and let its encoded latent be zr == De (XF; H.-1). At time step t, we obtain the counterfactual hidden state H. from the predicted counterfactual past (ŠC XCF pa,t-) The decoder recovers the true counterfactual at time step t almost surely:
XCF == क. (2,; AT = =X
Proof. We establish the result by induction over time steps f. Following the definition of the true counterfactual value X defined in (15) we denote the corresponding true hidden state H which updated using the true counterfactual values X and XCF pa.t-1-
Base case. Let to E I be the first interventional step after the context window. By the intervention rule,
XCF 7too and XCF "Ito
hence X XCF Moreover, since to is the first intervention, the histories up to to are factual both constructions, so
CF HCF,
66-1-
Induction step. Suppose that for time step until - 1, the estimated counterfactual history matches the true one, i.c.,
CF
C
Since the counterfactual hidden state a function of past history until t 1, this implies
#Cr
H
CF
(33)
where the left-hand side is updated recursively from the estimated counterfactual history, and the right-hand side denotes the true hidden state under the counterfactuals.
By the deterministic and invertible property of the flow for fixed conditioning input, we have
Φ ( Hi-1); H.-1) X
(34)
From the structural causal model, the factual sample is
x 5(Xt- Xpatt, U.),
17

---

Published as a conference paper at ICLR 2026
for some exogenous noise Up- Fixing the same Ut. and under the intervention do(X,- X Xpa.1 X pa.t ), the true counterfactual sample is
X²" =f f(x, x pa,t
(35)
By Proposition 4.3 the latent factual representation satisfies
z ox(: H[_1) - $ ( , H = 9(U)),
which depends only on U, and is invariant to (X, X p-,)-- and thus to H Consequently, under the same intervention do(X,- =x X p = XCF
pa,t we also have
Φ (Xf: B(-1) $ (f(( U2); Hz-1) = o(Ui) @b of x x U.;
(36)
and more simply, we have:
L (X.: 11m11) 2x ;111-1) ;
(37)
Combining (33 and (37.. we obtain
l :F:):: ( = (x, :
(38)
By 34 the right-hand side of (38) equals X,F,, Since the left-hand side of (38 is precisely the algorithm's counterfactual encoder-decoder procedure, follows that
x == 2 (2,; ( -
Therefore, the decoder recovers the true counterfactual at time step By induction, the claim holds for all t.
B CAUSAL SCM ASSUMPTIONS AND THEORETICAL COMPARISONS TO BGM
B. CAUSAL SCM ASSUMPTIONS
Let (xz) 7+1:T denote the potential trajectory under the intervention do(X, = 7z). where I C [K] x (I+1,..., T}. We assume the following standard conditions for the structural causal model on the given DAG:
Assumption B.1.
(I) Unconfoundedness. The exogenous noises (U,,) are mutually independent across (i, and all common causes of observed nodes are included in the DAG, i.e., no unobserved confounders for the interventions considered.
(2) Consistency. If the system is run under do(Xz == 72), then the observed trajectory equals the
corresponding potential trajectory under the regime actually applied: X x(1x)
(3) Positivity. For any (zi,e-, p in the support of X ). the conditional density
of X given these parents is positive in neighborhood of the intervention value %r
p (X., x X Zi,-- x Xpa(t)(- -ppa(),)-) 0.
(4) No Interference. For any two intervention schedules and J with the same assigned values on given index set, the corresponding potential trajectories coincide that set. Interventions affect other nodes only through the edges of the DAG.
Remark B.2. Assumption nodes are included in the DAG, so the exogenous noises (Ut,z) mutually independent across ti,t). This rules out time-varying latent confounders that simultancously influence multiple series. In settings where unobserved confounding is suspected, our identification and recovery guarantees for interventional and counterfactual distributions need not hold.
(1) imposes causal sufficiency: common causes of the observed
Extending DoFlow beyond the causal sufficiency assumption settings with latent confounders is an important direction for future work. Recent causal generative frameworks such as DeCaFlow (Almod) ovar et al. 2025) and Neural Causal Models (Xia et al. 2022) explicitly address hidden confounding
18

---

Published as a conference paper at ICLR 2026
problems. DeCaFlow augments causal normalizing flows with latent confounders and proxy variables to recover identifiable interventional and counterfactual queries. Besides, Neural Causal Models use neural SCMs with correlated exogenous variables and graphical constraints to decide and estimate identifiable counterfactuals from observational data. A promising direction is to combine such deconfounding mechanisms with our time-indexed CNF parameterization so that DoFlow can handle deconfounded temporal DAG.
B.2 THEORETICAL COMPARISONS TO BGM
Bijective Generation Mechanisms (BGM). The BGM framework (Nasr-Esfahany et al., 2023) shows that if the true structural mechanism f is in the BGM class (bijective/strielly monotone in the exogenous noise), then any learned mechanism f that (i) is also in the BGM class and (ii) matches the observed i.c.. for the same parents value X and action A,
f(Xi-, X.xt,,pu,t-a pp,tt U2) - 5(X,_, X Ui) i.e., (f(X,-, Xpa.- D) = (f(X,,, Xpatt D)
yields the same counterfactuals as Their result is -agee--gg,ottc, at the class level.
DoFlow. Our findings are specific particular methodology and architecture, this is methodol- ogy paper where we analyze continuous normalizing flows (CNFs). With the additional assumption (A3) specific to DoFlow, Proposition 4.3 shows that the DoFlow's 's encoded latent Z₁ is bijective in the proves that the encode- decode procedure recovers
exogenous noise Ug. Consequently, Corollary 4.5 the true counterfactual.
Relationships. Because of Proposition 4.3 (enabled by the CNF), DoFlow implements a bijective-in- noise mechanism and is thus comparable to BGM in the sense that f is bijective in U, However, they the same: BGM additionally requires observational distribution matching (as formalized in (A4)) to obtain class-level identifiability, whereas our proofs of Proposition 4.3 and Corollary 4.5 do assume (A4) and instead rely on the model-specific condition (A3).
Alternative route. The counterfactual recovery result (Corollary 4.5) can alternatively be established by imposing an additional assumption:
(A4) Observational matching. For each (X,-, Xpa,e-), the DoFlow-induced observational law matches the true one, i.e.
-f(x,
Under (A1), Proposition 4.3 establishes a continuously differentiable bijection g : U4 - Z (independent of (Xr-, Xpa.;-). ). so the induced mechanism f implemented by DoFlow is bijective in the exogenous noise, i.e., DoFlow lies in the BGM class. With the additional assumption (A4), we can therefore directly invoke the counterfactual recovery result of the BGM framework
al. 2023).
C
PRELIMINARIES ON CONTINUOUS NORMALIZING FLOWS
For better logical flow in the main text, we move the CNF preliminaries here. Continuous normalizing
flows (CNFs) have been successful for statistical sampling (Tian et 2024: Wu and Xie 2025) and counterfactual inference (Wu et al. 2025 Chao et al. 2024). This appendix section reviews CNF fundamentals and flow matchingg training for general data types.
C.1 NEURAL ODE AND CONTINUOUS NORMALIZING FLOW:
Neural ODE models the evolution of sample as the solution an ordinary differential equation (ODE). Concretely, in Rd, given an initial condition zo z00) at 0, the transformation to the output .i = x11) at s is governed by:
dr(s) ds
= v(z(s), s), s se [0,1],
(39)
19

---

Published as a conference paper at ICLR 2026
where Rd x [0.1] - IR is the velocity field parameterized by neural network. The time horizon is rescaled to 8 E [0, 1] without loss of generality.
Continuous Normalizing Flow (CNF) is class of normalizing flows in which the transformation of probability density is governed by a time-continuous Neural ODE. Let p(z, s) denote the marginal density of as)). Then p(ar, 8) evolves according the Liouville continuity equation implied by (39):
0,p(æ, V (p(#, 8) v(2, s)) 0, E 10,
(40)
where V. denotes the divergence operator.
When the Neural ODE is well-posed, it induces continuous and invertible map from the initial sample Zu to the terminal sample The inverse map is obtained by integrating (39 backward in time. This mechanism allows one to choose Iy ~ p(-, 0) as the data distribution and 21 p(: 1) as
a simple base (noise) distribution, typically [N, 7). For convenience, we write q(-) == p(., 1) for the base distribution. Throughout, we use the common choice q(-) = N(0,) I).
C.2 FLOW MATCHING
Flow Matching (FM) (Lipman et al. 2023) trains continuous normalizing flows without simulating trajectories by regressing the model velocity v(r(s), 8) toward a prescribed target field u(x()), The FM objective is
PPM E0.. [[lv((((), 8) (z((),
(41)
where u(x(s), s) is analytically specified.
Linear Interpolant: Because u(x; 8) is generally intractable unless we condition on the starting point xo. Conditional Flow Matching (CFM) (Lipman et al. 2023) was proposed to find a tractable velocity. Given an observed data zo and an endpoint æ1 (( from base distribution, we can choose an analytic interpolation between Zo and ø1, and define a reference path. An -rrrstrrion-rrreriing and simple choice the linear interpolant (Lipman et 2023):
(1 - s)z0 + (s min s))21,
(42)
where æ1 ~ N(0, 7) and amin is a small positive hyperparameter ensuring p(d, 0) ~ ( o T).
Setting min 0 recovers the strict linear path (1 s)ar + 881-
We treat á as a fixed reference trajectory that the learned flow is trained to track. Under 42 the associated reference velocity
dø
(1 Ominmaa 30-
(43)
ds
Training loss: CFM trains the flow v(((s),)) s) by directly regressing it onto the reference velocity field 43 The training objective is given by:
LCFM Be
dø
dØ
112
(44)
ds
D DATA SYNTHESIS AND METRICS DEFINITION
We define four types of causal DAG structures:
Tree, Diamond, Fully Connected Layer (FC-Layer), and Chain with skip connections (Chain). For each structure, we design both additive models and nonlinear, non-ad con-additive structural causal models.
In the simulations for different graph structures, the root node X1 (for Tree, Diamond, and Chain), and the root nodes X1, X2, X3 (for FC-Layer) are initialized over the interval [0, to] using a Chain process:
(45) Here, i E {1} for Tree, Diamond, and Chain, and 2,, 3) for FC-Layer. Besides, A. P, and Φ; denote the amplitude, period, and phase offset of the sinusoidal driver, and Uz, represents independent exogenous noise 10,
== = X1, + sin 2nd + +)) U. P
In the Tree, Diamond, and Chain structures, the hyperparameters for the root node X1,- namely A, P. and o. are set to 1.0, 20, and 25 respectively. For the FC-Layer, these are set to 1.5, 15, and 4x respectively.
20

---

Published as a conference paper at ICLR 2026
D.1 TREE
%
X
X
X
X X
Figure 4: Tree graph over 8 nodes. Exogenous variables UL. are omitted for clarity but exist for every node at each time Left: Full node-level causal structure between consecutive time, with all variables (X1.).... X8.1 present at each step. Right: Rolled-up (time-suppressed) view over different nodes {XI..... X Each arrow X, + X, (with i * j) denotes lag-1 temporal dependency X Xj.t that holds for all t. Both panels depict the same underlying structure.
We consider both additive and nonlinear, non-additive structural causal models:
Additive model:
X X 3,X Σ jEpa(i)
BingX + Un//4. (46)
Nonlinear, non-additive model:
X1,1 1.
B.X (|Ual + 0.5) +
Σ
Bie (47)
jEpa(i)
For the Tree graph, i I....,& (in order), and define edge-specific parent coefficients 3 only on the directed edges 1 2 (0.3), - (-0.3), 2 (0.3), 2 + (-0.3), 3 - 6 (0.5), 3 + (-0.5), and 7 (0.5), with 3ir- 0 for all other pairs.
the self-lag coefficients B E {0.5, 0.4, 0.2, 0.,
0.2, 0.2, for
D.2 DIAMOND
We consider both additive and nonlinear, non-additive structural causal models:
Additive model:
X1, =
- Σ
B + Un.
(48)
jEpa(1)
Nonlinear, non-additive model:
X = f(Xi,-1i Xpa6)))-1, Une) = exp(8xXx.011)
+
Σ
B
2 2+ Uke
jEpa(s)
(49)
21

---

Published as a conference paper at ICLR 2026
Ap
x
Figure 5: Diamond graph over 10 nodes. Exogenous variables U are omitted for clarity but exist for every node at each time t. Left: Full node-level causal structure between consecutive time, with variables 1X,,,..., X present at each step. Right: Rolled-up (time-suppressed) view over different nodes (X1,.... X10). Each arrow X; X, (with # j) denotes lag-1 1 temporal dependency X,, + Xj,t that holds for all t. Both panels depict the same underlying structure.
For the Diamond graph, we set the self-lag coefficients (B,,..., ;10) (0.5, 0.4, 0.3, 0.2, 0.1, 0.2, 0.2, 0.2, 0.2, 0.2), and define edge-specific parent coefficients Beey only on the directed edges >2 (0.3) 1-33 (0.3), 4 (0.3), - (0.3), 3->5 (0.5), 3 > (0.5), 4- (0.5), 8 (0.5), 5 - (0.5), 7- (0.5, 10 (0.5), and 10 (0.5), with Pie-j for all other pairs.
D.3 FULLY CONNECTED LAYER (FC-LAYER)
consider both additive and nonlinear, non-additive structural causal models:
Additive
X1, = f(Xie Xpatt-- - Σ
B + Uir-
(50)
jEpa(s)
Nonlinear, non-additive model:
Xi, = 1(Xix-1, Xpa(i), Un,e) BiXis-11 Σ
Biom X
jEpa(i)
? exp Σ
D -1 - tanh
Σ jEpa(t)
Bio Xzit-1) + Un
jEpa(í)
(51)
For the FC-Layer graph, we set 7 = 0.30 and the self-lag coefficients (3),..., 310) = (0.5, 0.6, 0.7, 0.4, 0.45, 0.5, 0.55, 0.3, 0.35. 0.4). We define edge-specific parent coefficients Bir densely between adjacent layers: for every hidden node i E (4,5, 6, 7) and input E 2,,2,3) we set 3ie_ 0.2 + 0.1 (j 1) (i.e., 0.2, 0.3, 0.4 for j 1,2,3), and for every output node i {8,9, 9, 10} and hidden E {4,5, 6, 7] we set Bie = 0.2, with Bic = 0 otherwise.
D.4 CHAIN LINEAR (CHAIN)
We consider both additive and nonlinear, non-additive structural causal models:
22

---

Published as a conference paper at ICLR 2026
X
XG
Figure 6: FC-Layer graph over 10 nodes. Exogenous variables U,, are omitted for clarity but exist for every node at each time l. Left: Full node-level causal structure between consecutive time, with all variables {Xi,t)..,, X10.e) present at each step. Right: Rolled-up (time-suppressed) view over different nodes ( dependency X + Xz, that holds for all t. Both panels depict the same underlying structure.
Each arrow X; + X, (with * j) denotes lag-1 temporal
time
Xooa
%
Figure 7: Chain graph over 50 nodes. Exogenous variables U4,; are omitted for clarity but exist for every node at each time Left: Full node-level causal structure between consecutive time, with all variables [X,,)..., Xso,t)
present at each step. Right: Rolled-up (time-suppressed) view Each arrow X; Xj (with & j) denotes lag-1 temporal dependency X + Xj that holds for all I. Both panels depict the same underlying structure.
different nodes (XI,..., X
Additive model:
X1, = (X Use)) = 3:X Σ jEpa(s)
B Uie
(52)
Nonlinear, non-additive model:
X1. = (.x,)),-X .U.))) = (3.54 1. 11-.11111111111111111 X1,1-1
(53)
23

---

Published as a conference paper at ICLR 2026
For the Chain structure, we set 3g = 0.6 for all i ,1...,,5,, 50, and draw independent Gaussian noises US ~ NN(0, a)) with 0; = 0.2. For each non-root node = 2.... 50, the immediate predecessor immediate-predecessoo coefficient is fixed as Bie-i-1 = 0.2, and the skip connection from the root is fixed as 3 0.7.
D.5 INTERVENTIONAL AND COUNTERFACTUAN SIMULATIONS
Each time series window has length T 120. We set the context window 90 and the forecasting window to T - 30. Interventions begin immediately after the context window (1.2....,
To rigorously evaluate DoFlow under challenging causal interventions, we design interventions that alter the root node(s) the causal DAG, since manipulating root causes propagates system-wide effects through all downstream variables. Each root node follows a cyclic structural equation (45).
To ensure that the intervention produces a fundamentally different dynamic pattern, construct the interventional trajectory by phase-shifting the root cycle by half a period, i.e., P/2. This yields:
Ni. = 9 + Asin P/2 + +0) + Usar tE
(54)
where i E {1) for Tree, Diamond, and Chain, and i E (1,2,3) for FC-Layer. This corresponds to complete inversion of the cyclic trend for the root nodes of each graph structure (including noise term create different intervention schedules across runs). Consequently, the interventional trajectory exhibits an antiphase pattern relative to the original, creating the most challenging regime for causal forecasting and counterfactual recovery.
Below, we illustrate how to simulate the true interventional and counterfactual trajectories using the SCMs defined above. For clarity, we present the procedure for the Chain, Diamond, and Tree structures, where the intervened root node is X1, The FC-Layer case follows identically, except that interventions are applied to root nodes E {1, 2, 3} as defined (54).
The interventional simulation proceeds by first sampling the exogenous noises {Uir+ I. R=2 The intervened values are then obtained by applying do(Xi,=-1:7 71.+1.7 and propagating forward with the sampled noises to generate the intervened trajectories {Xi++1:2)), 1K For the FC-Layer graph, the intervened parental includes { f...),,11,,23
Counterfactual simulation begins by recovering the exogenous noises (e(r+1122 from the observed factual future [XTr 1:T ],K using the structural causal models. Next, the intervention do(Xi. 7:.++1:7) is applied, and the system is propagated forward with the recovered noises to generate the counterfactual trajectories (XCF 1
6,7+1.77 1.T h- For the FC-Layer graph, the intervened
parental set includes {X 1x
D.6 METRICS
We evaluate model performance using both Root Mean Squared Error (RMSE) and Maximum Mean Discrepancy (MMD). Let the test batch size be B = 128, the ground-trute -ruth value for node at time t in batch b be xi. and the corresponding model prediction be -.
(b)
For each test sequence b. given the same context (x:(0) } generate N 100 realizations of both the model-estimated and the true observational//ntervaaa obtain more accurate evaluation metrics.
forecasting trajectories, in order
To ensure comparability across different scales time series, we apply standard normal scaling to
each batch over the forecasting window, using the mean μ and standard deviation (b) computed from its context window of node i. The generated time series dataset has a total length of 20,000, with stride of 1, resulting in 15,881 training samples (80%) and 3,881 testing samples (20%).
24

---

Published as a conference paper at ICLR 2026
D.6.1 RMSE
The RMSE of node for single realization is defined
T
B
RMSE;
Σ t=T+1 b=1
2 11)
(55)
B
For observational and interventional testing, we evaluate each method over 50 test batches of batch size B = 128. Within batch, RMSE is computed by (i) averaging over nodes = l,.... K, then (ii) averaging over N 100 stochastic realizations, and finally (ii) averaging over the 50 batches. The reported standard deviation (std) reflects variation across the 50 batches.
For counterfactual testing, each method produces a single deterministic counterfactual trajectory for given test instance, so N = 1. We evaluate performance by averaging RMSE over counterfactual test batches (each with B = 128), and report std across 50 batches.
D.6.2 MMD
To calculate MMD, we first flatten each trajectory (length - T, dimension D) into a vector in R (T)D, The sample size of both the true and the estimated trajectories (:7 JKK is BN, since for each batch b we simulate N = 100 realizations. The empirical MMD is then defined as
MMD
Σ
((Fa)Ta)+ + (BN)
k(Ts, Fv)
2 (BNY²
k(ra 3.),
(BN)(BN
0,0
aska
bydb'
(56)
where we use the Gaussian kernel
k(2, x²) exp
||--
2
and the bandwidth a is chosen via the pooled median heuristic.
D.6.3 CRPS
For probabilistic forecasting, we evaluate the quality of the predictive distribution using the Continu- ous Ranked Probability Score (CRPS) (Matheson and Winkler 1976). For a univariate predictive CDF F and a realized outcome E IR, the CRPS is defined as
CR = 1 (F(z) - 11f2 2 v3) dz,
(57)
where 1 1(-} } is the indicator function.
In our experiments, is represented by an ensemble of samples (a(h.-)) NN test sequence with realized value (b) We use the standard ensemble approximation:
for node i, time t, and
CRPS
(b) E
(b,)
(2)
(bim)
((h,n')
(58)
N
2N
report CRPS averaged over all test sequences b, time steps t E { and runs, analogous to the aggregation used for RMSE.
T}. nodes = 1,..., K,
Notably, CRPS is applicable only to probabilistic models. Therefore, we report CRPS results only for DoFlow, Dee, and MQF2.
E MODEL AND IMPLEMENTATION DETAILS
This section summarizes the main implementation and hyperparameter choices used all experi- ments.
25

---

Published as a conference paper at ICLR 2026
Common setup. All models are implemented in PyTorch, with continuous normalizing flows (CNFs) integrated using torchdiffedy odeint. Unless otherwise stated, we use: Adam opti- mizer with learning rate 10 batch size 128, and ReLU activations in the CNF networks. For each node i the DAG, we instantiate for our algorithm: (i) an LSTM-hased RNN encoder that takes as input the node value (and, when applicable, time features), and (ii) CNF whose velocity field is parameterized by 3-layer MLP with hidden width 64. All time series are normalized per window using a mean-variance scaler over the context window.
E. ADAPTING OBSERVATIONAL BASELINES FOR INTERVENTIONAL ROLLOUTS
For each baseline method, we enforce the causal factorization by training separate model for each node i, using input rolling context window consisting of the node history and its parents' histories, (X1,1=-1, X pa)(). 1:t- }, and predicting the next-step target X,, For architectures that natively support multi-horizon prediction, we use the same model class but set the prediction length to one step and apply recursively.
To predict interventional future under an intervention schedule I with values {)t,,), we perform multi-step rollout over the forecast window by iterating = + 1,.... 7 and forecasting nodes in topological order: if (i, I fix पi Ni. otherwise predict æi, from the baseline model using the latest rolled-out histories. The rolled-out predictions (including fixed intervention values) are then fed back as inputs for subsequent time steps, ensuring that interventions propagate to downstream nodes through the DAG.
E.2 SYNTHETIC SCM EXPERIMENTS
For the synthetic DAG experiments (Tree, Diamond, Chain, and FC-Layer graphs), we first generate long multivariate time series using the SCMs in Appendix D We then extract sliding windows for training and evaluation.
Each window has total length T 120, with context window of r 90 steps and forecasting horizon of -- = 30 steps. We construct windows with stride along the simulated sequence. The resulting windows split along the time index into 80% for training and 20% for testing, using the first 0.8N time stamps for training and the remaining 0.2.N for testing.
For all synthetic structures, we set up algorithm:
Per-node RNN encoder: 3-layer LSTM with hidden size 15.
Per-node CNF: MLP with hidden dimension 64 and 3 layers; the conditioning dimension is
cond_dim, 15 + |pa(1)) 15.
We train DoFlow for 15 epochs on the synthetic datasets. For interventional evaluation, we draw 100 flow samples per window to plot the confidence interval well to calculate the metrics.
E.3 HYDROPOWER FORECASTING
For the hydropower case study, we use processed Argonne hydropower data individual locations. From the raw time-stamped series, we extract windows of length
pontextllengte 180, prediction_lengiey = 30,
with stride 5 along time. Thus each training and test window contains 180 context time steps followed by a 30-step forecasting horizon.
We incorporate the time features, including the time of day and calendar features, to obtain time-feature dimension of TF We set up our algorithm:
Per-node RNN encoder: 3-layer LSTM with hidden size 15.
Per-node CNF: MLP with hidden dimension 64 and layers; the conditioning dimension is
cond_ din_ 15 |pa(f)] 15,
26

---

Published as a conference paper at ICLR 2026
including the node's ' own hidden state, and the hidden states of its parents.
For training (performed offline and loaded from checkpoints in the main hydropower experiments), use batch size 128, Adam with learning rate 10 10-3 and Nopochs = 15 with early stopping based on node-wise training loss.
E.4 CANCER TREATMENT EFFECT
For the cancer-treatment experiment, we use the public dataset from Bica et al. (220aa), preprocessed into patient-level sequences with actions and outcomes, We construct per per-patient windows of length
context len = 55, pred len
5,
where the context contains the past tumor volumes and treatments, and the prediction horizon contains the next 5 steps to be modeled under different treatment policies.
Here distinguish between outcome and action histories via separate RNNs in our algorithm:
Outcome RNN: 4-layer LSTTM with hidden size H 15 and input size 1 (normalized cancer volume).
Action RNN: 4-layer LSTM with hidden size H 15 and input size D 4 (chemo application, radio application, chemo dose, radio dose).
The CNF for the cancer outcome Y; is parameterized as a 3-layer MLP with hidden width 64, with conditioning dimension
cond_ dim 15 15,
corresponding to the concatenation of the current outcome hidden state and action hidden state.
We construct a training set patient windows via CancerTra IDatasse with batch size 128 and shuffle the patients at each epoch. We train the trnate-treatment model for 15 epochs, using the same per-window mcan-variance normalization as in the other experiments, All reported interventional and counterfactual metrics are computed by rolling out DoFlow over the 5-step prediction horizon under alternative (simulated) treatment policies, as detailed in Section F.4
F MORE EXPERIMENTAL RESULTS
F.1 COMPUTATIONAL COSTS
The number of training samples is 15,881 for each simulated datasets. We use a batch size of 128, and all experiments-iicluding training and sampling time comparisons-are conducted on single A100 GPU.
Table Comparison of model size.
Tree
Diamond
Layer
Chain
Hydropower
DoFlow
94, 664
121, 658
133, 946
646, 178
99, 592
GRU
36, 830
50, 106
732
124, 412
52, 026
TFT
860
100. 472
775,775
531, 026
95,318
TIDE
108.301
118. 332
118, 332
544, 539
122, 462
TSMixer
117,818
119.
118
550,
110. 118
DAR
73. 712
89, 280
89, 280
184, 280
84, 290
MQF2
128,677 677
160, 846
160, 846
702, 681
125, 459
27

---

Published as a conference paper at ICLR 2026
Table Comparison of training time per epoch, epochs convergence, and observational sampling time (for 1,000 forecast series over a 30-step horizon) across Tree, Diamond, and Layer structures.
Methods
Training Time / Epoch
Epochs to Conv.
Total Training Time
Sampling Time
DoFlow
42.3s
10
8.09min
10.25s
GRU
17.8s
15
5.03min
6.80s
TFT
21.0s
15
5.28min
11.07s
TiDE
29.3s
25
12n1
2.23s
TSMixer
27.9s
25
12. 1min
8,83s
DeepVAR
20.6s
20
11.65 11.65min
7.86s
MQF2
49.88
20
16.6min
19.61s
F.2 SYNTHETIC DATA EXPERIMENTS
Table 4: MMD for observational, interventional, and counterfactual time series forecasting across causal structures: Tree, Diamond, FC-Layer, and Chain. Reported values are averaged over 50 test batches, each containing 128 test series.
Tree
Dimmnnd
PC-layer
Chain
Addrine
NLNA
Addithe
NLNA
Addivv
NLNA
Addivv
MLNA
Ots.
Obs
Obs
Obs
3.00 "0.0P
018.a1
0.0.
0.00..
0.14.
0.17 14
0.01
1014
0.13..
101
010,,s
GRU
0.34
0.30
6.55.
0.0000
030100
0.1920
024
0.01
0.16
0..9..1
0.7..0
0.10..02
0.15
0.56.
0.31
TFT
m0.0..m
0.33a
6.20mm
0302
025
018
10200
022g
0.38gm
0.3
0,80
0,30
T0c
1000.
0.32mm
003em
10.22
0033gm
10.9.m
0.32
0022m
003em
018g/m
10.9.
0.35g
DapVAR
6.
0.00.s3
10.08
0.00.a
0.00u
0111
025..
008
023xa
1036..s
0,3
0.31.
0.1
0.34.
0.20,
M0F2
10.08.a:
0.30.
u.29.s
1000
036
0,95
030
0077.
0.00um
1021...
0.34.
0.17
0.22. 0.28
0.58
Table 5: CRPS for observational and interventional time-series forecasting across causal structures: Tree, Diamond, FC-Layer, and Chain. Reported values are averaged over 50 test batches, each containing 128 test series.
Tru
Diamand
FC-Layer
Chain
Adkiine
NENA
Addiine
NINA
Additiv
NENA
Addiin
NL.NA
Obs.
Obe
Obc
Obs.
In:
Obe
Obe.
020.
0.10.
0.7..
0010 0.4..0 on 0.24.
1038
Deep
0.37 1.
6,
0.0.0 e. u14g 0.9 (0.00 n 0.55 0.07g 0.28, M.
0.72
MQF2
10 10.19,
0.20g 0.542
6.10gm 1 0.77.m
033.m
03t.w 013em
0. 10322
0.29 0.38 0.32 0:46. 0.32
Table 6: RMSE for observational, interventional, and counterfactual time-series forecasting on the Chain causal structure.
Chain
Additive
NLNA
Obs.
Inl.
CF:
Obs.
Inl.
CF
DoFlow
0.61
0.66
0.30_
0.69
0.71
0.41::m 07
GRU
0.68
1.01
NA
0.80
1.21
NA
TFT
1.63
1.10
0.78
1.29
NA
TDE
0.65 m
1.07.17
0.7T..13
1.16
NA
TSMixer
0.67
1.09
NA 0,75
1.20
NA
VAR
0.62...2
0.97
NA 0.844.m
1.09
NA
MQF2
0.73
1.18
NA
0.90
1.30+
NA
28

---

Published as a conference paper at ICLR 2026
Nntevvndiional Forecasting
Axuneredi
M M
Figure 8: Left: "Tree" interventional forecasting results. Node X1,: is intervened. DoFlow provides 50% and 90% prediction intervals; the orange lines indicate true interventional future. Right: "Layer" counterfactual forecasting results. Nodes X1,ri X2,5, and X3,: are intervened. DoPlow provides a single forecast in green; orange lines indicate the true counterfactual future.
Counterfactual Porecasting
Figure 9: "Chain" (NLNA) counterfactual forecasting results. Node X1, is intervened. DoFlow
provides a single forecast in green; the orange lines indicate the true counterfactual future.
F.3 HYDROPOWER SYSTEM
Generator
Turbine
X
X
X
Transtomee
M&C
Figure 10: Hydropower system graph over 8 nodes. Exogenous variables UL. are omitted for clarity but exist for every node at each time Left: Full node-level causal structure between consecutive time, with all variables {X1, view over different nodes 1, Xs,3 Xx,
present each step. Right: Rolled-up (timessuprressed) Xs). Each arrow X; + X, (with &f j) denotes lag-I temporal dependency X; - X3, that holds for all t. Both panels depict the same underlying structure.
,X,,....,
29

---

Published as a conference paper at ICLR 2026
Reliable modeling and uncertainty quantification are increasingly central to modern signal processing (Dong et al., 2020; Nemani et al.] (2033; Gao et al. [023] Wu et al.) 2024a). Here, we highlight a new perspective enabled by DoFlow, causal reasoning in time series, through hydropower signal processing case study.
In Figure T0
1 and X2 denote the horizontal and vertical rotational vibrations of the hydraulic turbine, which directly drive the generator's ' operation. X3 represents an auxiliary turbine state that affected by X, and X2. The generator's horizontal and vertical dynamics are represented by X4 and Xs, while X6 captures the generator's 's current output delivered to the transformer (Xs) for voltage regulation and transmission the power grid. The Metering and Control (M&C) unit (X7) continuously monitors the turbine's 's performance and stability to ensure coordinated operation.
T the RMSE results for observational and interventional time-series forecasting in the hydropower system. For interventional forecast- ing, we use 12 real power outages where the root nodes X1 and X2 fail, causing the entire system to outage. In this setting, the root nodes are treated as intervened by the breakdown signals. The reported averages and standard deviations are computed over 12 runs with batch size B = 1 for the interventional case, in contrast to 50 runs with batch size = 128 for the observational case.
Hydropower System
Obs.
Int.
DoFlow
1.13
1.21
GRU
2.05_.32 32
2,45
TFT
1.82_.2
2,16+.41
TiDE
1.49_.24
2,08
TSMixer
1.51_22
2.11
DeepVAR
1.78_.26 26
2.39
Notably, the hydropower signals are highly un- stable, with turbine flow and generator readings often exhibiting abrupt jumps or burnouts without clear patterns. Consequently, all methods face diffi- culty in accurate prediction, and the relatively high RMSE values reflect this inherent challenge. Nev- ertheless, our model performs consistently better than others under these conditions.
MQF2
1.97_.24
2.62
Table 7: RMSE for observational and interven- tional time-series forecasting in the hydropower system.
Hydropower System Anomaly Detection
Hydropower System Anomaly Detection
Kir
Fir
rine
Mu
Figure 11: Anomaly detection by DoFlow on real power outages in the hydropower system (two segments shown).
For anomaly detection, we use DoFlow as a likelihood- scoring model over predicted future win- dows. Fix node and current time in the test sequence. We treat the observed past [X1...., X,) context, sample future window of length L. 30 minutes from DoFlow, {X X X. }.
and compute the conditional log-density that DoFlow assigns to this generated forecast trajectory (i.e., the model-implied confidence of its own prediction) using the change-of- -arriallss formula. Formally, this yields the node-wise anomaly score
e, (L) = log P0. (Sae+1:44 M..s)
where Li,+41:+ denotes the forecast window decoded from sampled base noise 21,411:t4 under DoFlow conditioned on H Here, the log- likelihooo is computed using the Eq. 13).. restricted to
30

---

Published as a conference paper at ICLR 2026
the window [++1,
(og f u=b+1
1+4 Σ
l:v. V. m. ( Hin-1) do (59)
log xeeeene (Bretnann m..)
with i -b (zi, H1)-1) and Hi the RNN summary of the past up to u- as defined in Section 2.3 Intuitively, (1) measures how typical the entire upcoming 30-minute trajectory of node i appears under the dynamics learned from normal-operation data by the model.
To obtain quantitative decision rule, first compute (t) all test windows that do not overlap any outages and collect the empirical distribution of these scores for each node i. Let v) denote the 10th percentile of this normal-score distribution. We then flag time locally abnormal for node whenever
(+(t) Fj-
To avoid reacting to isolated noise spikes, we declare an anomaly only if at least one node exhibits 5 consecutive minutes of low scores. i.c., there exists and a run of times f,ttl..... ++4 such that (T) 7h for all J in that run. The anomaly detection time for a given outage is the earliest at which this criterion is met; in Figure 11 the corresponding detection markers are plotted relative to the known outage time (minute 50 in each panel).
In Figure TT notably, DoFlow's 's log-probability output becomes abnormal well before the outage as early as 20 minutes prior (e.g., Xe,z in the right panel) and as late as 10 minutes prior. This allows anomalies to be detected in advance the true outage.
F.4 CANCER TREATMENT EFFECTS
We apply DoFlow to interventional forecasting of cancer tumor outcomes. For each test patient, the model observes the first 55 days of factual history {X,, Y,) ]t=1. | where X, {X.,, ) represents chemotherapy and radiotherapy assignments and dosages, and Y, denotes the tumor volume.
Notably, unlike our other experiments, here only the tumor-olumme node Y, is modeled with CNE The treatment actions X, (chemotherapy and radiotherapy assignments/dosages) are fully observed and directly provided every time step. We still use RNNs to summarize the joint history of past treatments X, and past tumor outcomes Y, into hidden state H,-1, which serves the conditioning variable for the CNF of Y,
During the forecasting window (days 36-62), the treatment schedule is replaced by one of ten pre- defined intervention plans T' {(i,t) X₃,, <- m where j indicates the j-th treatment plan and 21, ..., 10). At each time step DoFlow estimates the interventional tumor outcome by sampling from the learned flow model:
Y, ' He-1))
21 ~ ~(0
(60)
where 2b denotes the learned reverse flow conditioned on the recurrent hidden state 41-1, which encodes the patient's 's historical outcomes and past treatments. The hidden state is updated autoregres- sively using the newly generated Y, and the active treatments X, = m
The model implementations detailed in Appendix E.4
31

---

Published as a conference paper at ICLR 2026
Treatments
time
X
AE
X;
X
X
Figure 12: Cancer Treatment DAG over 8 nodes. Exogenous variables U are omitted for clarity but exist for every node at each time Left: Full node-level causal structure between consecutive time, with all treatment variables {X,+]] and cancer tumor outcome Y, present at each step. Right: Rolled-up (time-suppressed) view different nodes. Each arrow X; Y denotes lag-1 temporal dependency X + Y₁ that holds for all t. Both panels depict the same underlying structure.
Table 8: Normalized RMSE, for causal treatment effects cancer tumor outcome. At future step RMSE is computed across all patient-option pairs as. Column groups represent the chemotherapy and radiotherapy application budgets (Ner 7r) in various data simulation scenarios.
% 5,
5, %
0
To 0,
DoFlow
CRN
RMSN
MSM
DoFlow
CRN
DoFlow
CRN
RMSN
MSM
1.25%
2.43%
3.16%
6.75%
0.49%
1.08%
1.35%
3.68%
0.94%
I.54%
1.59%
3.23%
4
1.7%
2.83%
3.95%
7.65%
0.76%
1.21%
1.81%
3.84%
1.10%
1.81%
2.25%
3.52%
5
2.08%
3.18%
4.37%
7.95%
0.88%
1.33%
2.13%
3.91%
1.27%
2.03%
2.71%
3.63%
6 2.74%
3.51%
5.61%
8.19%
1.09%
1.42%
2.41%
3.97%
1.69%
2.23%
2.73%
3.71%
7
3.22%
3,93%
6.21%
8.52%
1.33%
1.53%
2.43%
4.04%
2.01%
2.43%
2.88%
3,79%
We compute the normalized root mean-squared error (RMSE) the 7-th step across patients and treatment options as
En=l (Ynu4r Y Em.tt7)
RMSE,
(61)
N ΣN Yastt
where N is the total number of patient-option pairs in the test set. Smaller NRMSE indicates more accurate estimation of causal treatment effects.
Table reports the normalized RMSE results for causal treatment effect estimation on cancer tumor outcomes. We compare our method with three established baselines: CRN (Counterfactual Recurrent Network) Bica et al.] (2020a), RMSN (Recurrent Marginal Structural Network) Lim (2018). and MSM (Marginal Structural Model) Mansournia al. (002)). The datasets follow construction in Bica et al. (2020a), where "Ye and 7r denote the treatment-applicatioo budgets for chemotherapy and radiotherapy, respectively. Some bascline results Table S are adopted directly from Bica et al. (2020a).
32


Reviews:
[{"review_id": "JrvEub88Jd", "reviewer": "Reviewer_sw48", "summary": "This paper proposes **DoFlow**, a time-conditioned continuous normalizing flow (CNF) defined over a causal DAG for multivariate time series. Each node’s flow is conditioned on recurrent states of the node and its parents, enabling an encode–decode procedure that supports observational, interventional, and counterfactual forecasting via abduct–action–predict. The training uses a conditional flow-matching loss; the model also exposes trajectory log-likelihoods used for automatic change-point/anomaly detection. On the theory side, under assumptions (independent noise, monotonicity in the noise, and a base-law condition on the encoded latent), the authors show that the encoded latent is a bijective function of the exogenous noise and derive a counterfactual recovery result for the encode–decode scheme; they relate this to BGM, noting their result does not require observational distribution matching. Empirically, they evaluate on synthetic DAGs (tree/diamond/chain/FC-layer) for interventional and counterfactual rollouts, and on a real hydropower system, reporting interventional/observational RMSEs and demonstrating early outage detection from log-densities (often 10–20 minutes in advance).", "questions": "* Following weakness 2, do you have any results for the calibration of the output intervals?\n* Did you consider to replacing the RNN by a small transformer or just an attention layer? It would be more costly at prediction time, but would be interesting to see if it improves performance, and there could maybe even have some speed-ups in the case of a large past time series with a small forecast, due to parallelization.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "* Encoder–decoder CNF over a causal DAG neatly instantiates **abduction–action–prediction**, so observational, interventional, and counterfactual forecasts all come from one coherent mechanism.\n* **Impressive empirical results**: strong observational/interventional accuracy and compelling real-world demos (incl. early anomaly/outage detection), plus robustness across diverse synthetic DAGs.\n* **Theory that matches the use case**: clear assumptions lead to bijectivity-in-noise and a **counterfactual recovery** result that directly supports the encoder–decoder design.\n* **Clear, well-structured presentation** that makes the model easy to follow and reproduce.", "weaknesses": "* **Positioning vs prior causal–flow work could be clearer.** The paper does a good job situating the CNF side (Neural ODEs, flow matching, encode–decode) but is lighter on contrasting with prior *causal* modeling using flows (e.g., works like *Causal Normalizing Flows: From Theory to Practice*, Javaloy et al.). A short subsection disentangling what comes from causal literature, what from flows, and what’s **novel here** would help readers map contributions more precisely.  \n* **Prediction intervals lack calibration evidence.** Interventional figures show 50%/90% bands, but I couldn’t find quantitative **coverage** or calibration metrics (e.g., empirical coverage, interval width/ACE, Winkler score, CRPS) to support the claim. Consider reporting coverage on synthetic data (where ground truth is known).\n* **Discussion section is a bit too sparse in my opinion** Even though the results seem good, I think a small discussion on particular on the limitations of this work would be interesting.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "pIayABdJsA", "reviewer": "Reviewer_dDVY", "summary": "The objective of this work, as expressed in Equation (9), is to obtain the counterfactual forecast $X_{\\tau+1:T}^{\\mathrm{CF}}$, given the interventions $\\mathrm{do}(X_I := \\gamma_I)$ applied to the predictor $X_{1:\\tau}$ and the corresponding factual forecast $X_{\\tau+1:T}^{\\mathrm{F}}$, along with optional conditional variabels.  The authors trained a mapping that connects the current state to its past states and parent variables, based on a directed acyclic graph (DAG), while optionally incorporating conditional variables.  Using these factors as inputs, the authors employed flow matching to generate the current state as a probability distribution over possible outcomes. This process is autoregressive. To support the training of the entire pipeline, the authors used simulated data generated from four causal patterns, encompassing both linear and nonlinear additive models. After the entire framework is trained, given the essential inputs and specified interventions,  it can generate the corresponding counterfactual forecast distribution.", "questions": "1. Please address the presentation issues mentioned above.\n\n2. According to the mechanism of generative AI, particularly the flow-matching approach adopted in this work, the output represents only a sample from a potential distribution. Would this introduce additional uncertainty that complicates interpretability? Even if the analysis relies solely on the probability density function (PDF) to quantify significance, there remains a potential limitation regarding practical applicability: during training, small perturbations in past states may lead to extreme events in the future, and vice versa. Since such cases are statistical outliers, how can the authors justify that the trained flow-matching model is capable of capturing these rare patterns? If not, there is a risk that the forecasted extremes may be erroneously dismissed when relying exclusively on the PDF. \n\n3. This work appears to focus on utilizing past information, possibly generated by applying interventions to an existing dataset. However, what is the advantage of introducing an additional corresponding past state (perhaps the original one before intervention)? Why not treat this simply as a standard forecasting problem? Although the authors claim that causality plays a role, this assertion is not sufficiently justified. Moreover, since one past state could result from different interventions applied to various factual cases, would incorporating such additional causal structures introduce greater complexity or instability compared to a conventional forecasting framework?\n\n4. In Appendix B.5, the details on how the interventions were constructed are still missing. To demonstrate the practical value of this work, especially in real-world scenarios involving complex patterns, it is essential to justify that both the model design and the simulated data can effectively support training for complex pattern recognition. Unfortunately, the current version of the paper does not provide sufficient evidence to substantiate this key point.\n\n5. The reproducibility of this work is another major concern. The paper lacks essential records of hyperparameter settings, which is particularly important given the inherent outcome uncertainty of the generative AI model employed. Furthermore, the details of how the simulation data were constructed and how the model was trained to adequately capture complex real-world patterns, so as to perform as claimed in the case study, are not clearly documented. These omissions make the reproducibility and reliability of the reported results questionable.\n\nI would consider raising my scores if the above concerns are adequately addressed.", "limitations": "", "rating": 2, "confidence": 5, "soundness": 2, "presentation": 1, "contribution": 2, "strengths": "1. This paper introduces an unexplored research direction, as claimed by the authors, which could be of potential value. \n\n2. The technical pipeline presented in this work appears to be complete. \n\n3. Producing the probability density function (PDF) instead of a single counterfactual forecast is an insightful approach.", "weaknesses": "1. Several key concepts are not clearly introduced at the beginning. For example, in the Introduction section, it is unclear what is meant by “interventional question.” How are “binary or discrete, fixed-time actions, treatment focus …” defined? What does “extrapolating a counterfactual path and contrasting it” mean? Most importantly, the key concept of “counterfactual” should be clearly illustrated from the beginning. What does “these works do not encode causal DAGs or simulate system-wide trajectories under interventions” mean?\n\n2. In the Introduction, I would prefer the authors to directly explain why the research problem is important, rather than simply stating that it has not been studied and is important in many fields without sufficient evidence or support. The authors are also advised to briefly outline the key challenges of the problem and how this paper addresses them, instead of letting the related work take up a large portion of the Introduction.\n\n3. Why are the contributions placed in the Related Work section instead of the Introduction? Why does the Related Work section repeat content that has already been mentioned in the Introduction? In addition, the authors have not conducted a sufficiently thorough review of related work. Does the proposed method fall under an existing category of related work? How does it differ from similar studies?\n\n4. It is recommended that the authors introduce the concept of a causal DAG in the Preliminary section. It should explain more on how a multivariate time series can be represented as a directed graph. It sounds like paths in a graph, but can we still name them as time series?\n\n5. For the technical representation, the authors did not explain the meaning of the notation $\\gamma$ in equation (8), and did not explain the difference between notations X and x. The meaning of h should be explained early for equation (10), where it first appeared. Equation (11) is not clear. Is it a concatenation of vectors of scales or a mapping? What is the specific use of the conditions? In equation (11), I suggest the authors to remind the readers the meaning of s clearly: the time step of flow matching or the time series. The authors also did not explain the meaning of K nodes. How do we construct the nodes? What are they representing?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "dZ7amIcGQG", "reviewer": "Reviewer_jniY", "summary": "This paper proposes _DoFlow_ a framework based in NeuralODEs and Recurrent Neural Networks (RNN), to estimate not only observational forecasting distributions, but also interventional distributions and counterfactuals, opening the door to reliable decision making in decision making with time-series.\n\nDoFlow models each variable in a known causal graph with its own NeuralODE, and conditions generation on a recurrent state that summarizes the past of that variable and its parents. This lets the model similate future trajectories in topological order. This allows to forecast observational series, interventional distributions and individual counterfactual outcomes. Because the model is a continuous normalizing flow, it also assigns likelihoods to predicted trajectories, which they use for anomaly and change-point detection. The authors provide a recovery result guaranteeing correct counterfactual reconstruction under certain assumptions, and they report experiments on synthetic causal time series and a real hydropower system. DoFlow matches or improves standard forecasting error compared to strong baselines, and uniquely supports interventional simulation, counterfactual rollout for each individual sequence, and early anomaly warnings before system failures.", "questions": "- I do not really understand the implication that _counterfactual explanations_ in this work. Although the term _counterfactual_ is a common term in causal literature, counterfactual explanations are non-causal concepts that look for perturbations in the input to get the desired output. Is it necessary to include that in the introduction? I think those papers only introduce noise to the reader.\n\n- Related with the previous point, I have the same question with causal discovery literature. Although causal discovery is a much more related topic, why do you include that in the introduction, if you assume a known causal graph? I think the introduction is scarce, in general very noisy, and the objective of the work is not clear.\n\n- I am surprised because in the related work the authors name some models that solve interventional and counterfactual queries, but they do not name two models that are the state-of-the-art in this field and are based in normalizing flows, thus being very connected with the technologies of this paper: [1, 2]. Is there any reason for not including them?\n\n- In general, some related works are missing or they lack discussion. For example, haven't authors considered to include [3]? Or can authors discuss the limitations of [4]? Those are only some examples, but more discussion have to be added in the introduction/ related work.\n\n- Should the paper of continuous normalizing flows not be cited in line 115?\n\n- Although section 4 is well written in general, $\\Phi_\\theta$ should be defined before introducing them in the algorithms. Otherwise, the reading is more difficult. Instead, it is defined in section 5 (`line 322`).\n\n- Assumption A3 requires $Z\\sim N(0,1)$, but is that distribution necessary? Is not enough with having Z independent from $H_{t-1}$, independently of the base distribution of $Z$?\n\n- Everything that the DAG is assumed to be true, there is a question that arises: what if the DAG is not well specified? I think that some discussion about this fact would help.\n\n- Have the authors thought of exploring the violation of 'all confounders observed'? Some work on this have been explored in causal generative models: [5, 6]\n\n> References\n\n[1] Javaloy, A., Sánchez-Martín, P., & Valera, I. (2023). Causal normalizing flows: from theory to practice. Advances in Neural Information Processing Systems, 36, 58833-58864.\n\n[2] Khemakhem, I., Monti, R., Leech, R., & Hyvarinen, A. (2021, March). Causal autoregressive flows. In International conference on artificial intelligence and statistics (pp. 3520-3528). PMLR.\n\n[3] Liu, Y., Sun, Y., & Lim, J. H. (2023, June). Counterfactual Dynamics Forecasting–a New Setting of Quantitative Reasoning. In Proceedings of the AAAI Conference on Artificial Intelligence (Vol. 37, No. 2, pp. 1764-1771).\n\n\n[4] Shenghao Wu, Wenbin Zhou, Minshuo Chen, and Shixiang Zhu. Counterfactual generative models\nfor time-varying treatments. In Proceedings of the 30th ACM SIGKDD Conference on Knowledge\nDiscovery and Data Mining, pages 3402–3413, 2024\n\n[5] Almodóvar, A., Javaloy, A., Parras, J., Zazo, S., & Valera, I. (2025). DeCaFlow: A Deconfounding Causal Generative Model. arXiv preprint arXiv:2503.15114.\n\n[6] Kevin Muyuan Xia, Yushu Pan, and Elias Bareinboim. Neural Causal Models for Counterfactual Identification and Estimation. In The Eleventh International Conference on Learning Representations, ICLR 2023", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- Preliminaries about NeuralODE  and Flow matching are clear. Additionally, the method section: 'Time conditioned flow on a causal DAG'  is very well written and the method is easy to understand. The method is based on state of the art deep learning techniques, and all the proposed strategy is sound and practical. I find the model useful for future reseach. However, I would like to see a complete scheme of the implementation of the framework: how many NeuralODEs are needed, how the RNN is trained and combined with the Flow, etc.\n\n- The paper presents theoretical guarantees about the recoverability of counterfactuals, which are sound and mostly (see weaknesses) clear, although I have not checked the proofs of the Appendix in detail.\n\n- The results on synthetic data support the theoretical insights in interventional forecasting and in counterfactual estimation. Also, experiments on real-world data demonstrate an impresive performance, not only in estimating interventions, but also in estimating anomaly detection via likelihoods. Experiments are systematic, extensive and complex enough.", "weaknesses": "- Neither the problem or the solution proposed are well specified in the introduction. There are many unconected paragraphs that do not provide a connected story line about the motivation, the **limitations of existing work** or the strategy that this paper proposes. For example, `line 96`says: ''motivated by these gaps...'', which gaps? Also, more information about the interventions that the model is suposed to handle has to be provided. What does '...counterfactual queries across the entire system'(`line 97-98`) mean?\n\n- I have a minor concern about the theory section, although I probably miss something. First of all, I do not see a condition of invertibility between X and U. Assumption A2 requires monotone $f(\\cdot, U_t)$. However, that is invertible only if the function is continuous and also if X is continous. How can we recover the value of the exogenous variable if the function is not bijective? I think that bijectivity should be required explicitly.\n\n- There are implicit assumptions that have not being discussed in the paper:\n1) Unconfoundedness / causal sufficiency. I.e., that all confounders are observed. If there exist time-varying confounders that are not observed, I do not see how the model can recover interventional distributions and counterfactuals.\n2) Positivity. If the observational data is not rich enough, that is, if the intervention is done in a point in which the density of the observational distribution is not zero, the model would extrapolate the causal effects.\n3) no interference and consistency: the classical SUTVA assumptions.\n4) Regularity conditions: I think that the functions of the SCM have to meet some regularity conditions in order this theory to hold.\n\n\n> Summary\n\nAlthough I have some concerns---specially with the introduction and the related work section---and many questions, I find the paper a good contribution to causal inference in time series. I think the employment of neuralODEs to causal time series forecasting, with the inclusion of DAGs, which allows for interventional and counterfactual inference, is great step forward in causal forecasting. For all the previous reasons, I will recommend to accept the paper, and I will be happy to raise my score if my concerns are covered.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "XUZ5ALCRWC", "reviewer": "Reviewer_Pjb9", "summary": "This paper proposes a novel architecture based on Normalizing Flows for time-series forecasting. The model is designed to handle three types of forecasting queries: observational (standard prediction), interventional (predicting the effect of an action), and counterfactual (predicting what would have happened under a different action, given a specific observed outcome).\nWhile the idea is interesting and the authors provide an appropriate theoretical analysis of the method, I have some concerns about specifically the presentation, the empirical evaluation, and the motivation for counterfactual simulations that I would like to be addressed before recommending the paper for acceptance.", "questions": "- As noted in the specified weaknesses. \n\n- What happens if the causal graph is not fully known or partly wrongly specified? Is DoFlow able to recover from such misspecifications?\n\n- How would you position your method against works such as [1] ?\n\n- Could you imagine under which conditions (or data distributions) DoFlow is struggling? A discussion on this could help the reader assess the applicability of your method to their problem. \n\n¹ https://arxiv.org/abs/2402.09891", "limitations": "", "rating": 4, "confidence": 2, "soundness": 2, "presentation": 2, "contribution": 3, "strengths": "1.  **Novelty:** The core idea of using a unified Normalizing Flow-based framework to tackle observational, interventional, and counterfactual time-series forecasting is interesting and novel.\n2.  **Theoretical Grounding:** The paper provides a theoretical analysis of the proposed method, which adds rigor to the contribution.", "weaknesses": "1. **Motivation and Practical Utility of Counterfactual Forecasting:** My most fundamental concern is the practical relevance of the counterfactual simulations as defined in the paper. The method hinges on conditioning on the exogenous noise `z` from an observed trajectory to ask \"what if?\" questions. However, the motivation for why a practitioner or decision-maker would need to perform a simulation conditioned on this specific noise instance is unclear.\n    - Could the authors provide a concrete, practical decision-making scenario where this type of counterfactual simulation provides unique, actionable information that cannot be obtained from simpler interventional forecasts?\n    - How does this concept translate to models that are not probabilistic (e.g., deterministic models)? Is intervening on a variable `X` sufficient for counterfactual reasoning if there are no further noise components?\n\n\n2. **Insufficient Experimental Details and Reproducibility:** The empirical evaluation is missing critical information, making it difficult to assess the validity of the results.\n    - **Hyperparameters and Model Specifications:** There are no details on the hyperparameters for the proposed model or for any of the baselines. Information on architecture, training duration, learning rates, etc., is essential. E.g. 908-910 is extremely vague and should be further specified.\n    - **Dataset Details:** The paper lacks information on the amount of data used for training, validation, and testing for each experiment, nor does it provide information on whether such a split was actually used. For the real-world power outage data, it is unclear how the model was trained, as there is presumably only one historical trajectory. Was the power outage event included in the training data, or was the time series split? This needs further clarification.\n    - **Intervention Details:** The nature of the interventions is not specified. Were the interventions in-distribution or out-of-distribution? How many variables were affected simultaneously? This context is crucial, especially given that the baselines appear to perform strikingly worse than the proposed method, even when trained on only the correct Parent set. A detailed description of the experimental setup is needed to understand this performance gap.\n\n3. **Unclear Explanation of Anomaly Detection Results:** The application to preemptive anomaly detection is highlighted as a key result, but is not sufficiently explained.\n    - How exactly is the detection performed? What quantity is measured, and what threshold defines an \"abnormal\" reading?\n    - The claim of *preemptive* detection is extraordinary. Why can the model detect an anomaly before it occurs? Is this consistent with the physical reality of the system, or is it an artifact of the experimental setup (e.g., information leakage)? A more detailed explanation of the underlying mechanism is necessary to properly assess this result.\n\n4. **Clarity and Connection of Paper Sections:** Certain parts of the paper seem disconnected or lack support.\n    - **L60 (Causal Discovery):** The discussion of the causal discovery literature feels vague and disconnected from the paper's core contribution, as the method fundamentally assumes a known and fixed causal graph (DAG). The authors should either better integrate this discussion or remove it.\n    - **Section 4.4:** This section appears unsupported and feels out of place, as it is not explained how exactly the anomaly detection can be performed.\n    - **Methodological Limitations:** It should be explicitly noted and discussed that the proposed model does not account for instantaneous causal effects within a single time step, which is a significant assumption in many real-world systems.\n    - **Table 1:** The table is overflowing and quite small. Please consider updating the formatting", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)