# Structured Learning of
Compositional Sequential Interventions

Jialin Yu
University College London
University of London
&Andreas Koukorinis
University College London
University of London
Nicolò Colombo
Royal Holloway
University of London
&Yuchen Zhu
University College London
University of London
&Ricardo Silva
University College London
University of London

###### Abstract

We consider sequential treatment regimes where each unit is exposed to
combinations of interventions over time. When interventions are described by
qualitative labels, such as “close schools for a month due to a pandemic” or
“promote this podcast to this user during this week”, it is unclear which
appropriate structural assumptions allow us to generalize behavioral
predictions to previously unseen combinatorial sequences. Standard black-box
approaches mapping sequences of categorical variables to outputs are
applicable, but they rely on poorly understood assumptions on how reliable
generalization can be obtained, and may underperform under sparse sequences,
temporal variability, and large action spaces. To approach that, we pose an
explicit model for _composition_ , that is, how the effect of sequential
interventions can be isolated into modules, clarifying which data conditions
allow for the identification of their combined effect at different units and
time steps. We show the identification properties of our compositional model,
inspired by advances in causal matrix factorization methods but focusing on
predictive models for novel compositions of interventions instead of matrix
completion tasks and causal effect estimation. We compare our approach to
flexible but generic black-box models to illustrate how structure aids
prediction in sparse data conditions.

##  1 Contribution

$D^{1}_{n}$$D^{2}_{n}$$D^{3}_{n}$$D^{4}_{n}$$\beta_{n}$$X^{1}_{n}$$X^{2}_{n}$$X^{3}_{n}$$X^{4}_{n}$

Figure 1: Within unit $n$, interventions $D_{n}^{t}$ interact with (latent)
random effect parameters $\beta_{n}$ to produce behaviour $X_{n}^{t}$
represented as a dense graphical model with square vertices denoting
interventions (Dawid, 2021). Further assumptions will be required for the
identifiability of the impact of interventions and their combination,
including temporal impact.

Many causal inference questions involve a treatment sequence that varies over
time. In the (discrete) time-varying scenario, an _unit_ $n$ is exposed to a
sequence $D_{n}^{1},D_{n}^{2},D_{n}^{3},\dots$ of _treatments_ controlled by
an external agent (Hernán and J, 2020, Ch. 19). Each $D_{n}^{t}$ can be
interpreted as a cause of the subsequent behavior
$X_{n}^{t},X_{n}^{t+1},\dots$, and its modeling amenable to the tools of
causal inference (e.g., Robins, 1986, 1997; Murphy, 2003; Didelez et al.,
2006; Chakraborty and Moodie, 2013; Zhang and Bareinboim, 2020). For instance,
each $D_{n}^{t}$ can take values in the space of particular drug dosages which
can be given to an in-patient, or combinations of items to be promoted to a
user of a recommender system. Here, patients and users are the units, and
$X_{n}^{t}$ denotes a measure of their state or behavior under time $t$.
Dependencies can result in a dense causal graph, as in Fig. 1. Further
background on this classical problem is provided in Appendix A.

#### Scope.

We introduce _an approach to model the effects of a combinatorial and sparse
space of intervention sequences, including the bursts of non-stationarity that
the application of an intervention brings in._ We consider _the case where
each $D_{n}^{t}$ behaves as if it was randomized_ to focus on predictive
guarantees on large, possibly underspecified, sequential action spaces. We
will consider neither experimental design nor reinforcement learning, instead
focusing on _the problem of behavioral forecasting under hypothetical future
interventions given past controlled interventions_.

#### Challenges and prior work.

Much of the sequential intervention literature in AI and statistics considers
generic models111We include here models which explicit parameterize causal
contrasts in longitudinal interventions, such as the blip models of Robins
(1997). to map sequences $D_{n}^{1:t}:=[D_{n}^{1},D_{n}^{2},\dots,D_{n}^{t}]$
to behavior $X_{n}^{t}$ and beyond, sometimes exploiting strong Markovian
assumptions or short sequences and small action spaces (Chakraborty and
Moodie, 2013; Sutton and Barto, 2018). An unstructured black-box for
sequences, such as recurrent neural networks and their variants (e.g.
Hochreiter and Schmidhuber, 1997; Cho et al., 2014; Vaswani et al., 2017), can
learn to map a sequence to a prediction. This is less suitable when
interventions are categorical and sparse, with most entries corresponding to
some baseline category of “no changes”. That is, there is no obvious choice of
smoothing or interpolation criteria to generalize from seen to unseen strings
$D_{n}^{1:T}$, with the practical assumption being that we see enough
variability to reliably perform this interpolation. There is an increased
awareness of problems posed by large discrete spaces among pre-treatment
covariates (Yao et al., 2019; Zeng et al., 2024), which is of a different
nature as losses can be minimized with respect to an existing natural
distribution, as opposed to the interventional problem where the design
variables are allowed to be chosen and so sparsity conditions can easily
result in very atypical test cases of interest. Several papers address large
intervention spaces (mostly for non-sequential problems) without particular
concerns about non-smooth, combinatorial identifiability (Kaddour et al.,
2021; Nabi et al., 2022; Saito et al., 2023, 2024). A few more recent papers
address explicit assumptions of identifiability directly in combinatorial
categorical spaces by sparse regression (Agarwal et al., 2023a) or energy
functions Bravo-Hermsdorff et al. (2023), but without a longitudinal
component. More classically, tools like the do-calculus (Pearl, 2009) and
their extensions (Correa and Bareinboim, 2020) can be used with directed
acyclic graphs (DAGs) to infer combinatorial effects, although they often
restrict each intervention variable to target a single random variable (e.g.
Aglietti et al. (2020)). With strong assumptions and computation, this may
include latent variables Zhang et al. (2023). Alternative representation
learning ideas for carefully extrapolating to unseen interventions are
presented by (Saengkyongam et al., 2024) in the non-sequential setting.

In this paper, we consider the conservative problem of _identification
guarantees_ for the effect of sequential combinations of interventions which
may never have co-occurred jointly in training data, _formalizing assumptions_
that allows for the transfer and recombination of information learned across
units. These are not of the same nature as identifiability of causal effects
from observational data (Pearl, 2009) or causal discovery (Spirtes et al.,
2000), but of _causal extrapolation_ : even in the randomized (or
unconfounded) regime, it is not a given that the distribution of a particular
potential outcome $X(d^{1:t})$ can be reliably inferred from limited past
combinations of interventions without an explicit structure on the causal
generative model. We will consider the case where interventions have non-
stationary effects: once unit $n$ is exposed to $D_{n}^{t}=d_{n}^{t}$, then
$X_{n}^{t+\Delta}$ will have dynamics that can be affected not only by the
choice of $d_{n}^{t}$ but also the difference $\Delta$, resulting in dense
dependency structures such as the one in Figure 1. As multiple treatments are
applied to the same unit $n$, there will be a _composition_ of effects that
will depend not only on the labels of the interventions, but also on the order
by which they take place. If a substantive number of choices of interventions
at time $t$ is available, a large dataset may have no two sequences
$D_{n}^{1:T}$ and $D_{n^{\prime}}^{1:T}$ taking the same values, even when
considering only the unique treatment values applied $n$ and $n^{\prime}$
regardless of ordering and time stamps.

#### Problem statement.

Let each individual unit $n$ be described by time-invariant features $Z_{n}$
and a pair of time-series $(\\{X_{n}^{t}\\},\\{D_{n}^{t}\\})$. Each
$D_{n}^{t}$ is a categorical variable representing an _intervention_ (also
called _treatment_ , or _exposure_), encoded by values in $\mathbb{N}$, and
sequentially unconfounded with the system by randomization or assumption.
$D_{n}^{t}$ describes an intervention which causes changes to the distribution
of $X_{n}^{t:\infty}$, the time-series of state variables starting at time
$t$. The special value $D_{n}^{t}=0$ denotes an “idle” or “default” treatment
that can be interpreted as “no intervention”. Given a training set of units
$(Z_{n},\\{X_{n}^{t}\\},\\{D_{n}^{t}\\})$ and observations up to a time point
$T$, the goal is to predict future sequence of states for an unit $n^{\star}$
under a hypothetical intervention $D_{n^{\star}}$. We denote it as
$X_{n^{\star}}^{T+1:T+\Delta}(d_{n^{\star}}^{1:T+\Delta})$, where
$V^{t_{1}:t_{2}}$ denotes a slice of a discrete-time time series of random
vectors $V$ from $t_{1}$ to $t_{2}$, inclusive; $V(d^{1:t})$ is the potential
outcome (Imbens and Rubin, 2015; Hernán and J, 2020; Pearl, 2009) of variable
$V$ under interventions $D^{1}=d^{1},D^{2}=d^{2},...,D^{t}=d^{t}$. This
process takes place under stable unit treatment value assumption (SUTVA)
(Imbens and Rubin, 2015), meaning that interventions $D_{n}^{t}$ only affect
unit $n$.

Each intervention $D_{n}^{t}$ is allowed to have an impact on the whole future
series $t,t+1,t+2,\dots$. It can correspond to an instantaneous shock (“give
five dollars of credit at time $t$”) or an action that takes place over an
extended period of time (“promote this podcast from time $t$ to time $t+5$”):
we just assume the meaning to be directly encoded within arbitrary category
labels $0,1,2,3,\dots$. In a real-world scenario where a same action can be
applied more than once, different symbols in $\mathbb{N}$ should be used for
each instance (e.g., 1 arbitrarily standing for “give five dollars”, and 2 for
“give five dollars a second time” and so on).

In Section 2, we provide an account of what we mean by compositionality of
interventions, with explicit assumptions for identifiability. In Section 3, we
describe an algorithm for likelihood-based learning, along with approaches for
predictive uncertainty quantification. Further related work is covered in
Section 4 before we perform experiments in Section 5, highlighting the
shortcomings of black-box alternatives.

##  2 A Structural Approach for Intervention Composition

For simplicity, we will assume scalar behavioral measurements $X_{n}^{t}$, as
the multivariate case readily follows. Consider the following conditional mean
model for the potential outcome $X_{n}^{t}(d_{n}^{1:t})$,

$\displaystyle\mathbb{E}[X_{n}^{t}(d_{n}^{1:t})~{}|~{}x_{n}^{1:t-1}(d_{n}^{1:t-1}),z_{n}]=(\phi_{n}^{t})^{\mathsf{T}}(\beta_{n}\odot\psi_{n}^{t})=\sum_{l=1}^{r}\phi_{nl}^{t}\beta_{nl}\psi_{nl}^{t},$
(1)

where $\odot$ is the elementwise product,
$X_{n}^{t}(d_{n}^{1:t^{\prime}})=X_{n}^{t}(d_{n}^{1:t})$ for any
$t^{\prime}\geq t$ (future interventions do not affect past outcomes), and
$Z_{n}(d_{n}^{1:t})=Z_{n}$ for all $t$. We do not explicitly condition on
$D_{n}^{1:t}$ in most of our notation, adopting a convention where potential
outcome indices $d_{n}^{1:t}$ are always in agreement with the (implicit)
corresponding observed $D_{n}^{1:t}$. We define
$\phi_{nl}^{t}:=\phi_{l}(x_{n}^{1:t-1}(d_{n}^{1:t-1}),z_{n})$ as evaluations
of basis functions $\phi_{l}(\cdot)$, which for now we will assume to be
known. Each individual $n$ has their individual-level coefficients
$\beta_{n}$. Causal impact, as attributed to $D_{n}^{1:t}=d_{n}^{1:t}$, is
given by

$\psi_{nl}^{t}:=\prod_{t^{\prime}=1}^{t}\psi_{l}(d_{n}^{t^{\prime}},t^{\prime},t),$
(2)

where function $\psi_{l}(d,t^{\prime},t)$ captures a trajectory in time for
the effect of $D_{n}^{t^{\prime}}=d$, for $t^{\prime}<t$. We consider two
variants. The first is a _time-bounded_ variant which, for $d>0$, defines
$\psi_{l}(d,t^{\prime},t):=w_{dl,t-t^{\prime}}$, where $w_{dl,t-t^{\prime}}$
is a free parameter for $0\leq t-t^{\prime}<k_{d}$, otherwise
$w_{dl,t-t^{\prime}}=w_{dl,k_{d}-1}$. The intuition is that the effect of
intervention level $d$ has no shape constraints but it cannot affect the past
and it must settle to a constant after a chosen time window hyperparameter
$k_{d}>0$. The default level $d=0$ has no effect and no free parameters, with
$\psi_{l}(0,t^{\prime},t)=1$. The second variant is motivated by real-world
phenomena where causal impacts diminish their influence over time and result
in a new equilibrium (Brodersen et al., 2015). We consider a _time-unbounded_
shape with exponential decay, parameterized as

$\psi_{l}(d,t^{\prime},t):=\sigma(w_{1dl})^{t-t^{\prime}}\times
w_{2dl}+w_{3dl},$ (3)

where $\sigma(\cdot)$ is the sigmoid function, so that $0<\sigma(w_{1d})<1$.
As $t$ grows, $\sigma(w_{1d})^{t-t^{\prime}}$ goes to 0. $w_{3d}$ is the
stationary contribution of an intervention at level $d$. In particular, for
$d=0$, we define $w_{10l}=w_{20l}=0$ and $w_{30l}=1$, while for $d>0$ we have
that $w_{1dl}$, $w_{2dl}$ and $w_{3dl}$ are free parameters to be learned. By
defining $k_{d}$ as the _dimensionality of intervention level $d$_, we adopt
the convention $k_{d}=3$ for the prescribed time-unbounded model of Eq. (3).

###  2.1 Representation Power and Compositional Warping

Eq. (1) is reminiscent of tensor factorization and its uses in causal modeling
(Athey et al., 2021; Agarwal et al., 2023b), where the primary motivation
comes from the imputation of missing potential outcomes taking place on a pre-
determined period in the past. From a predictive perspective, it can be
motivated from known results in functional analysis (such as Proposition 1 of
(Kaddour et al., 2021)) that allows us to represent a function $f(x)$ by first
partitioning its input into $(x_{a},x_{b})$ and controlling the approximation
error given by an inner product of vector-valued functions
$f(x_{a},x_{b})\approx f_{a}^{\mathsf{T}}(x_{a})f_{b}(x_{b})$, in practice
choosing the dimensionality of the vectors by a data-driven approach. Although
those results apply typically to smooth functions with continuous inputs (such
as the Taylor series approximation, which relies on vectors of monomials with
coefficients constructed from derivatives), discrete inputs such as e.g.
$D_{n}^{1:T}$ at a fixed dimensionality $T$ can be separated from the other
inputs as
$f(n,d^{1:T},x^{1:T},z)=f_{d}^{\mathsf{T}}(d^{1:T})f_{nxz}(n,x^{1:T},z)$.
Here, similarly to classical ANOVA models, $f_{d}(\cdot)$ is a one-hot
encoding vector for the entire (exponentially sized) space of possible
$d_{n}^{1:T}$ trajectories, and $f_{nxz}(\cdot,\cdot,\cdot)$ returns a vector
of corresponding length.

As we do not want to restrict ourselves to fixed $T$ (and $N$), not to mention
an exponential cost of representation, we will resort to fixed-dimensional
embeddings of sequences and identities. In particular, we design
$f_{d}(d^{1:t})$ to have a finite output dimensionality $r_{g}$. Inspired by
sequential hidden representation models such as recurrent neural networks, we
define $f_{dl}(d^{1:t}):=g_{l}(d^{1:t},t,t)$, for $l\in 1,2,\dots,r_{g}$, via
the following recursive definition on $0\leq t^{\prime}\leq t$:

$g_{l}(d^{1:t^{\prime}},t^{\prime},t):=h_{l}(g_{l}(d^{1:t^{\prime}-1},t^{\prime}-1,t),d^{t^{\prime}},t^{\prime},t),$
(4)

such that $h_{l}(v,0,t^{\prime},t)\equiv v$ for any $(v,t^{\prime},t)$ (that
is, intervention level $D_{n}^{t^{\prime}}=0$ does nothing) and
$h_{l}(v,d,0,t)\equiv 1$ for any $(d,v,t)$ (base case). A representational
choice
$h_{l}(v,d,t^{\prime},t):=h_{l_{1}}(v)^{\mathsf{T}}h_{l_{2}}(d,t^{\prime},t)$
requires a choice of dimensionality $r_{h_{l}}$ for the inner vectors
$h_{l_{1}},h_{l_{2}}$. Recursively applying (4) in this product form can still
be done relatively efficiently, but choosing $r_{h_{l}}=1$ for all $l$
simplifies matters as we apply again the inner product formula to
$f_{d}^{\mathsf{T}}(d^{1:T})f_{nxz}(n,x^{1:T},z)$ to get to a factorization
form analogous to Eq. (1). More formally:

###### Proposition 1

Let $f(n,d^{1:T},x^{1:T},z):=f_{d}^{\mathsf{T}}(d^{1:T})f_{nxz}(n,x^{1:T},z)$,
where function sequences $f_{d}(d^{1:T})$ and $f_{nxz}(n,x^{1:T},z)$ are
defined for all $T\in\mathbb{N}^{+}$ and have codomain $\mathbb{R}^{r_{g}}$.
Furthermore, define $f_{dl}(d^{1:T})$ as in Eq. (4) with $r_{h}=1$, and the
$l$-th entry of $f_{nxz}(n,x^{1:T},z)$ as
$f_{nxzl}(n,x^{1:T},z):=u_{l}(n)^{\mathsf{T}}v_{l}(x^{1:T},z)$, each of the
latter having codomain $\mathbb{R}^{r_{s^{\prime}}}$ and where
$v_{l}(\cdot,\cdot)$ is defined recursively for any $T$, similarly
$g_{l}(\cdot,\cdot,\cdot)$. Then there exists some integer $r$ and three
functions $a,b$ and $c$ with codomain $\mathbb{R}^{r}$ so that
$f(n,d^{1:T},x^{1:T},z)=\sum_{l=1}^{r}a_{l}(n)\times b_{l}(x^{1:T},z)\times
c_{l}(d^{1:T},T)$. Moreover $c_{l}(d^{1:T},T)=\prod_{t=1}^{T}m_{l}(d^{t},t,T)$
for functions $m_{l}:\mathbb{N}^{3}\rightarrow\mathbb{R}$.

Proofs of this and next results are given in Appendix B. The above shows that
if we start with an assumption that the ground truth function takes some
desired form, then there is some finite $r_{g}$ for which our combination of
basis functions can parameterize the ground truth function for arbitrary $T$,
analogous to in Eq. (1). On the other hand, we show a companion result below,
that if we hold $T$ fixed, there is some fixed value for $r_{g}$, which
depends polynomially on $T$, such that any measurable ground truth function
can be approximated by our combination of some appropriate choice of basis
functions. This is reminiscent of the line of work on matrix factorization
(e.g., Appendix A of Agarwal et al. (2023c)).

###### Proposition 2

Given a fixed value of $T$, suppose we have a measurable function
$f:\bigcup_{t=2}^{T}\mathcal{X}^{1:t-1}\times\bigcup_{t=1}^{T}\mathcal{D}^{t}\times\mathcal{Z}\rightarrow\mathbb{R}^{d}$,
$t\leq T$, where $\mathcal{X}^{1:t-1}$ and $\mathcal{Z}$ are compact,
$\mathcal{D}^{t}$ is finite. For $d=1,\cdots,d_{max}$, let
$m_{d}:\bigcup_{t=1}^{T}\mathcal{D}^{t}\rightarrow\\{1,\cdots,T\\}$ with
$m_{d}(\mathbf{d})=t\text{ if }\exists t\;s.t.\;\mathbf{d}_{t}=d\;\text{else
}0$. For every $T$ there exist $r_{g}=O(T^{d_{max}+1})$ such that for any
$\epsilon$, there are measurable functions
$\phi_{l}:\mathcal{X}^{1:t-1}\times\mathcal{Z}\rightarrow\mathbb{R}^{d}$ and
$\psi_{ld^{\prime}}:\mathcal{D}^{t}\rightarrow\mathbb{R}$ and real numbers
$\beta_{l}$ for $l=1,\cdots,r_{g}$ and $d=1,\cdots,d_{max}$,
$\hat{f}(\mathbf{x},\mathbf{d},z)=\sum_{l=1}^{r_{g}}\phi_{l}(\mathbf{x},z)\beta_{l}\prod_{d=1}^{d_{max}}\psi_{ld^{\prime}}(m_{d^{\prime}}(\mathbf{d}))$
such that

$\displaystyle\text{sup}_{\mathbf{x}\in\bigcup_{t=2}^{T}\mathcal{X}^{1:t-1},\mathbf{d}\in\bigcup_{t=1}^{T}\mathcal{D}^{t},z\in\mathcal{Z}}\left|f(\mathbf{x},\mathbf{d},z)-\hat{f}(\mathbf{x},\mathbf{d},z)\right|\leq\epsilon.$
(5)

#### Notes.

Eq. (1) and Eq. (2) are motivated by the success of approaches for
representation learning of sequences, which assumes that we can compile all
the necessary information from the past using a current finite representation
at any time step $t$. However, knowledge of Propositions 1 and 2 further imply
that, for
$f(n,d^{1:T},x^{1:T},z):=\mathbb{E}[X_{n}^{t}(d_{n}^{1:t})~{}|~{}x_{n}(d_{n}^{1:t-1}),z_{n}]$,
the representation of functions in a multilinear combination of adaptive
factors (Eq. (1)) is motivated by more fundamental results in functional
analysis. Our assumption that no value of $d$ (other than the baseline 0) can
be used more than once has mostly an empirical motivation, as this would imply
unrealistic assumptions about the same intervention having the same effect
when applied multiple times. Ultimately, we assume that the factors
$h_{l_{2}}(d,t^{\prime},t)$ leading to Eq. (2) are time-translation invariant,
being a function of $d$ and $t-t^{\prime}$ only. While it is possible to
generalize beyond translation invariant models, as commonly done in the
synthetic controls literature (Abadie, 2021) that inspired the causal matrix
factorization ideas (Athey et al., 2021; Agarwal et al., 2023b, c), this would
complicate matters further by requiring a fourth factor into the summation
term of Eq. (1) to encode absolute time. Models for time-series forecasting
with parameter drifts can be tapped into and integrated with our model family,
but we leave these out as future work.

#### Interpretation.

At the starting of the process, where all units are assumed to be at their
“natural” state (assuming $D_{n}^{1}=0$), we can interpret
$\eta_{n}^{1}:=\beta_{n}$ as a finite basis representation, with respect to
$\phi_{n}^{1}$, of
$\mathbb{E}[X_{n}^{1}(0)~{}|~{}Z_{n}]=(\phi_{n}^{1})^{\mathsf{T}}\eta_{n}^{1}$.
Each $\eta_{nl}^{1}$ can be seen as a latent feature of unit $n$. Effect
vector $\psi_{n}^{t}$ defines a _warping_ of $\eta_{n}^{t-1}$ that describes
how treatments affect the behaviour of an individual with the resulting
$\eta_{n}^{t}:=\eta_{n}^{t-1}\odot\psi_{n}^{t}$. The modifier $\psi_{n}^{t}$
can be interpreted as reverting, suppressing or promoting particular latent
features that are assumed to encode all information necessary to reconstruct
the conditional expectation of $X_{n}^{t}(d^{1:t})$ from the chosen basis
(notice that intervention level 0 implies $\eta_{n}^{t}=\eta_{n}^{t-1}$ since
we define $\psi_{l}(0,t^{\prime},t)=1$). This presents a sequential warping
view of _compositionality_ of interventions.

###  2.2 Identification and Data Assumptions

Assume for now that functions $\phi$ are known. As common in the literature of
matrix factorization, assume also that we have access to some moments of the
distribution. In particular, for given $N$ data points and $T$ time points, we
have access to
$\mathbb{E}[X_{n}^{t}(d_{n}^{1:t})~{}|~{}x_{n}^{t-1}(d_{n}^{1:t-1}),z_{n}]$.
We will impose conditions on the realizations of $X_{n}^{1:T}$ and values
chosen for $d_{n}^{1:T}$, as well as values of $N$ and $T$ as a function of
$r$, in order to identify each $\beta_{n}$ and the parameters of
$\psi_{l}(d,t^{\prime},t)$ for all intervention levels $d$ of interest. The
theoretical results presented, which describe how parameters are identifiable
from _population_ expectations and _known function spaces_ determined by a
prescribed basis $\phi$ of _known and finite dimensionality_ $r$, will then
provide a basis for a practical learning algorithm in the sequel.

###### Assumption 1

For a given individual $n$, we assume there is an initial period
$T_{0_{n}}\geq r$ with “no interventions” i.e., $D_{n}^{t}=0$ for all $1\leq
t\leq T_{0_{n}}$.

###### Assumption 2

let $A_{n}$ be a $(T_{0_{n}}-1)\times r$ matrix where each row
$t=1,2,\dots,T_{0_{n}}-1$ is given by
$\phi_{n}(x_{n}^{1:t}(0^{1:t}),z_{n})^{\mathsf{T}}$. We assume that $A_{n}$ is
full rank with left pseudo-inverse $A_{n}^{+}$.

###### Proposition 3

Let $b_{n}$ be a $(T_{0_{n}}-1)\times 1$ vector with entries
$\mathbb{E}[X_{n}^{t+1}(d_{n}^{1:{t+1}})~{}|~{}x_{n}^{1:t}(d_{n}^{1:t}),z_{n}]$.
Under Assumptions 1 and 2, $\beta_{n}$ is identifiable from $A_{n}$ and
$b_{n}$.

###### Assumption 3

For a given intervention level $d$, assume $\psi_{l}(d,t^{\prime},t)\neq 0$
for all $l,t,t^{\prime}$ and that there is at least one time index $t$, and
one set $\mathcal{N}_{d}$ containing all units $n$ where $d_{n}^{t+1}=d$, such
that: (i) $N_{d}:=|\mathcal{N}_{d}|\geq r$; (ii) $\forall
n\in\mathcal{N}_{d},d_{n}^{t+2}=\dots=d_{n}^{t+k_{d}-1}=0$, where $k_{d}$ is
the dimensionality of the intervention level $d$.

###### Assumption 4

Let $n_{1},n_{2},\dots,n_{N_{d}}$ index the elements of a set of units
$\mathcal{N}_{d}$. For $t^{\prime}=1,2,\dots,k_{d}-1$, let $A_{dt^{\prime}}$
be a $N_{d}\times r$ matrix where each row $i=1,\dots,N_{d}$ is given by
$\phi_{n_{i}}(x_{n_{i}}^{1:t+t^{\prime}-1}(d_{n_{i}}^{1:t+t^{\prime}-1}),z_{n_{i}})^{\mathsf{T}}$.
We assume that each $A_{dt^{\prime}}$ is full rank with left pseudo-inverse
$A_{dt^{\prime}}^{+}$.

This leads to the main result of this section:

###### Theorem 1

Assume we have a dataset of $N$ units and $T$ time points
$(d_{n}^{1:T},x_{n}^{1:T},z_{n})$ generated by a model partially specified by
Eq. (1). Assume also knowledge of the conditional expectations
$\mathbb{E}[X_{n}^{t}(d_{n}^{1:t})~{}|~{}x_{n}^{t-1}(d_{n}^{1:t-1}),z_{n}]$
and basis functions $\phi_{l}(\cdot,\cdot)$ for all $l=1,2,\dots,r$ and all
$1\leq n\leq N$ and $1\leq t\leq T$. Then, under Assumptions 1-4 applied to
all individuals and all intervention levels $d$ that appear in our dataset, we
have that all $\beta_{n}$ and all $\psi_{l}(d,\cdot,\cdot)$ are identifiable.

Notice that nothing above requires prior knowledge of all intervention levels
which will exist, and could be applied on a rolling basis as new interventions
are invented. There is no need for all units to start synchronously at $t=1$:
the framing of the theory assumes so to simplify presentation, but the only
requirement is that each unit is given a “burn-in” period of at least $T_{0}$
steps and that each new intervention level $d$ is applied to at least $r$
units which have not been perturbed recently by relatively novel interventions
(that is, units $\mathcal{N}_{d}$ which, if they were given a past
intervention not given to anyone else before, then that perturbation should
have taken place enough steps in the past so that all recent but past $\psi$
functions have been identified).

##  3 Algorithm and Statistical Inference

This section introduces a learning algorithm, and well as ways of quantifying
uncertainties in prediction. We also allow for the learning of adaptive basis
functions $\phi$. As Eq. (1) does not define a generative model, which will be
necessary for multiple-steps-ahead prediction, for the remainder of this
section we will assume the likelihood implied by
$X_{n}^{t}(d_{n}^{1:t})=f(x_{n}^{1:t-1}(d_{n}^{1:t-1}),z_{n})+\epsilon_{n}^{t}$.
Here, the conditional mean $f(\cdot,\cdot)$ is given by Eq. (1) and
$\epsilon_{n}^{t}\sim\mathcal{N}(0,\sigma^{2})$, an error unaffected by
$D_{n}$ and other variables. Hence, $\sigma^{2}$ requires no further needs for
identification results. In general, if our likelihood is defined via a finite
set of estimating equations, we can parameterize it in the multilinear form
analogous to (1) and repeat the analysis of the previous section. In practice,
the functionals (such as the conditional expectations of $X_{t}$) used in an
estimating equation are unknown. Matrix factorization methods can be applied
directly by first performing a smoothing method on the data to get plug-in
estimates of these functionals Squires et al. (2022), but we will adopt
instead a likelihood-based approach.

###  3.1 Algorithm: CSI-VAE

We treat each $\beta_{n}$ as a random effects vector, giving each entry
$\beta_{nl}$ an independent zero-mean Gaussian prior with variance
$\sigma^{2}_{\beta}$. Along with all hyperparameters, we optimize the
(marginal) log-likelihood by gradient-based optimization. With a Gaussian
likelihood, the posterior and filtering distribution of each $\beta_{n}$ can
be computed in closed form, but in our implementation we used a black-box
amortised variational inference framework (Kingma and Welling, 2013; Mnih and
Gregor, 2014; Rezende et al., 2014) anyway that can be readily put together
without specialized formulas and is easily adaptable to other likelihoods. We
use a mean-field Gaussian approximation with posterior mean and variances
produced by a gated recurrent unit (GRU) model (Chung et al., 2014) composed
with a multilayer perceptron (MLP). The approximate posterior parameters at
sample with $t$ time steps are therefore
$\mu_{q,\beta,n}:=\text{MLP}(\text{GRU}_{\eta_{\beta,1}}(d_{n}^{1:t},x_{n}^{1:t},z_{n})),$
and
$\log\sigma_{q,\beta,n}:=\text{MLP}(\text{GRU}_{\eta_{\beta,2}}(d_{n}^{1:t},x_{n}^{1:t},z_{n}))$.
Prediction for $X_{n}^{t+1:t+\Delta}$ is done by sampling $M$ trajectories,
where for each trajectory we first sample a new $\beta_{n}$ from the mean-
field Gaussian approximation. We set each forward $\hat{X}_{n}^{t+i}$ to be
the corresponding marginal Monte Carlo average. Each basis vector
$\phi_{n}^{t}$ is parameterized via another GRU model such that
$\phi_{n}^{t}:=\text{MLP}(\text{GRU}_{\phi}(x_{n}^{1:t-1},z_{n}))$. Although
the theory in the previous section relies on known basis functions, we could
train this GRU up to an threshold time point $T_{0}$ and freeze it from that
point, each $\beta_{n}$ being defined as the unit coefficient vector for this
basis. In practice, we found that doing end-to-end learning provides a modest
improvement, and this will be the preferred approach in the experiments. We do
find that sometimes it is more stable to only use time series up to time point
$T_{0}$ for the amortized $\beta$ inference. Hence, we adopt this pipeline for
any results reported in this paper, unless otherwise specified. We call our
method _Compositional Sequential Intervention Variational Autoencoder (CSI-
VAE)_.

###  3.2 Distribution-free Uncertainty Quantification

Conformal prediction (CP) Vovk et al. (2005); Fontana et al. (2023) provides
prediction intervals with coverage guarantees. The intervals are computed
using a calibration set of labeled samples and include the future samples with
non-asymptotic lower-bounded probability. We consider Split Conformal
Prediction in two setups.

  1. 1.

Hold-out predictions. We have a set of _historical users_ , $n=1,\dots,N$,
whose behaviour has been observed until time $t+\Delta$. The task is to
predict the behaviour at time $t+\Delta$ of a _new user_ , $n=N+1$, who has
been observed up to time $t$, i.e. to predict $X_{N+1}^{t+\Delta}$ given
$X_{N+1}^{1:t}$ and $D_{N+1}^{1:t}$. If we assume we have used the history of
the new user $X_{n}^{1:t}$, $n=1,\dots,N+1$, to train the model, calibration
and test samples are exchangeable.

  2. 2.

Next-intervention predictions. We have observed the behaviour of _all users_ ,
$n=1,\dots,N+1$, up to time $t$ and aim to predict the effects of the next
intervention, which happens at time $t+1$ for all users, i.e. $D_{n}^{t+1}\neq
0$, holding $D_{n}^{t+2:t+\Delta}=0$. The task is to predict
$X_{n}^{t+\Delta}(d_{n}^{1:t+\Delta})$ for $n=1,\dots,N+1$, but in what
follows we will drop the potential outcome notation to keep notation lighter,
referring to observed $D$. Calibration and test are _not exchangeable_ ,
because i) the joint distribution after time $t$,
$P_{T}\sim(X_{n}^{t+\Delta},X_{n}^{1:t},D_{n}^{1:t})$ may be different from
the one before time $t$,
$P_{C}\sim(X_{n}^{t^{\prime}+\Delta},X_{n}^{1:t^{\prime}},D_{n}^{1:t^{\prime}})$,
$t^{\prime}<t-\Delta$ and ii) we only used
$X_{n}^{t^{\prime}},D_{n}^{t^{\prime}}$, $t^{\prime}=1,\dots t$, for training.

CP algorithms are applied on top of given point-prediction models. In our
case, the underlying model is $f$, which predicts the expected user behaviour
at time $t+\Delta$ given the user history up to time $t$, i.e.
$\hat{X}_{n}^{t+\Delta}:=f(D_{n}^{1:t},X_{n}^{1:t},Z_{n})\approx{\mathbb{E}}[X_{n}^{t+\Delta}~{}|~{}D_{n}^{1:t},X_{n}^{1:t},Z_{n}]$.

#### Setup 1.

Calibration and test scores,
$\\{S_{n}=|X_{n}^{t+\Delta}-f(D_{n}^{1:t},X_{n}^{1:t},Z_{n})|\\}_{n=1}^{N}$
and $S_{N+1}=|X_{N+1}^{t+\Delta}-f(D_{N+1}^{1:t},X_{N+1}^{1:t},Z_{n})|$ are
exchangeable, i.e. ${\rm Prob}(S_{1},\dots,S_{N}+1)={\rm
Prob}(S_{\sigma(1)},\dots,S_{\sigma(N+1)})$ where $\sigma$ is any permutation
of $\\{1,\dots,N+1\\}$. The Quantile Lemma, e.g. Lemma 1 of Tibshirani et al.
(2019), implies the prediction interval

$\displaystyle C=\\{x\in{\cal
X},|x-\hat{X}^{t+\Delta}_{N+1}|\leq\hat{Q}_{\alpha}\\}$ (6)
$\displaystyle\hat{Q}_{\alpha}=\inf_{q}\left\\{\sum_{n=1}^{N}{\bf
1}(|X^{t+\Delta}_{n}-\hat{X}^{t+\Delta}|\leq q)\geq n_{\alpha}\right\\},\quad
n_{\alpha}=\lceil(1+N)(1-\alpha)\rceil$

is _valid_ in the sense it obeys

$\displaystyle{\rm Prob}\left(X_{N+1}^{t+\Delta}\in C\right)\geq 1-\alpha,$
(7)

where the probability is over the calibration and test samples.

#### Setup 2.

The prediction intervals defined in (6) may not be valid, i.e. (7) may not
hold, because the calibration and test samples,
$S_{n}^{t^{\prime}}=|X_{n}^{t^{\prime}+\Delta}-f(D_{n}^{1:t^{\prime}},X_{n}^{1:t^{\prime}},Z_{n})|$
with $t^{\prime}\leq t$ and $t^{\prime}>t$, $n=1,\dots,N$, are not
exchangeable. Theorem 2 provides a bound on the coverage gap, i.e. a measure
of violation of (7), under the assumption that the distribution shift is
controlled by a perturbation parameter, $\epsilon>0$.

###### Theorem 2

Assume we have $N$ calibration samples,

$\displaystyle
S_{n}^{t^{\prime}}=|X^{t^{\prime}}_{n}-f(D^{1:t_{n}}_{n},X^{1:t_{n}}_{n},Z_{n})|,\quad
t^{\prime}=t_{n}+\Delta<t,\quad n=1,\dots,N$ (8)

where $t_{n}$ is the time user $n$ experienced the last intervention before
$t$. Assume there exists $\epsilon>0$ such that, for all $n$,

$\displaystyle
p_{T}(S_{n}^{t+\Delta})=(1-\epsilon)p_{C}(S_{n}^{t+\Delta})+\epsilon
p_{\delta}(S_{n}^{t+\Delta}),$ (9)

where $p_{T}$ and $p_{C}$ are the (unknown) densities of the test and
calibration distributions, $p_{\delta}$ is a bounded arbitrary shift density,
and $p_{min}=\min_{n=1,\dots,N}p_{C}(S_{n}^{t^{\prime}})>0$. Then,

$\displaystyle{\rm Prob}\left(X_{N+1}^{t+\Delta}\in C\right)\geq
1-\alpha-\frac{1}{p_{min}}\frac{\epsilon}{1-\epsilon}.$ (10)

In the proof, we use the likelihood-ratio-weighting approach of Tibshirani et
al. (2019) to obtain the empirical test distribution from the calibration
samples and bound its quantile from below. The statement follows from standard
inequalities on the mean and variance of the $\epsilon$-perturbed
distribution. We prefer this approach to conformal prediction adaptive models
for time series, e.g. Gibbs and Candes (2021) or Angelopoulos et al. (2024),
for two reasons: i) adaptive schemes have asymptotic coverage guarantees that
can not be used to estimate the uncertainty on a single time step and ii)
optimised density estimates are a byproduct of our prediction model.

##  4 Further Related Work

Brodersen et al. (2015) leverages Bayesian structural time-series models to
estimate causal effects, and motivates our exponential decay model in Eq. (3).
Unlike Brodersen et al. (2015), who focus on single interventions, our model
explicitly addresses the complexity arising from sequential interventions,
providing a more detailed perspective on the dynamic interplay of treatments
over time. Several approaches for conformal prediction in causal inference
have emerged in recent years (e.g., Tibshirani et al., 2019; Lei and Candès,
2021), including matrix completion Gui et al. (2023) and synthetic controls
Chernozhukov et al. (2021). Our focus has not been on individual nor average
treatment effects, but directly on expected potential outcomes, similarly to
the work on causal matrix completion Athey et al. (2021); Agarwal et al.
(2023b, c); Squires et al. (2022). Unlike the matrix completion literature, we
focused on prediction problems, out-of-sample for both units and time steps.
Moreover, the matrix completion literature is usually framed in terms of
_marginal_ expectations $\mathbb{E}[X_{n}^{t}(d)]$, as opposed to conditional
expectations (notice that Agarwal et al. (2023c) considers covariates, but
those are included as part of the generative model). Marginal models have some
advantages when modeling multiple-step-ahead effects Evans and Didelez (2024),
but they also involve complex computational considerations that we leave for
future work. Finally, there is a rich literature on confounding adjustment
where exogeneity of $D_{n}^{t}$ cannot be assumed, see Hernán and J (2020) for
a textbook treatment. Some of the causal matrix completion literature takes
the perspectiven that factors such as $\beta_{n}$ could be interpreted as
latent factors conditional on which ignorability could be assumed. However, we
consider this assumption hard to justify due to the unknown nature of
$\beta_{n}$. We can rely on standard approaches of sequential ignorability
Robins (1997) to justify our method in the absence of randomization.

##  5 Experiments

We run a number of synthetic and semi-synthetic experiments to evaluate the
performance of the CSI-VAE approach. In this section, we summarize our
experimental set-up and main results. The code for reproducing all results and
figures is available online222The code will be released upon acceptance, and
is available in the supplementary material.; in Appendix C, we provide a
detailed description of the datasets and models; in Appendix D, we present
further analysis and more results; and in Appendix E we present a demo for
uncertainty quantification results.

Table 1: Main experimental results, averaged mean squared root error over five
different seeds.

| Full Synthetic | Semi-Synthetic Spotify
---|---|---
Model | T+1 | T+2 | T+3 | T+4 | T+5 | T+1 | T+2 | T+3 | T+4 | T+5
CSI-VAE-1 | $36.53$ | $41.46$ | $41.73$ | $41.12$ | $41.32$ | $68.23$ | $82.94$ | $83.53$ | $81.97$ | $79.63$
CSI-VAE-2 | $97.80$ | $118.25$ | $117.79$ | $127.25$ | $135.03$ | $253.85$ | $312.53$ | $305.08$ | $303.68$ | $302.83$
CSI-VAE-3 | $138.78$ | $164.02$ | $141.71$ | $132.59$ | $125.55$ | $757.94$ | $937.07$ | $800.55$ | $704.66$ | $634.72$
GRU-0 | $229.72$ | $269.66$ | $220.95$ | $208.30$ | $188.43$ | $215.42$ | $260.65$ | $193.41$ | $137.20$ | $117.06$
GRU-1 | $230.76$ | $270.83$ | $220.93$ | $208.33$ | $184.92$ | $223.61$ | $269.69$ | $205.91$ | $141.53$ | $126.36$
GRU-2 | $93.73$ | $101.03$ | $118.01$ | $88.53$ | $132.28$ | $154.18$ | $187.42$ | $177.96$ | $133.36$ | $127.58$

|
---|---
|
Figure 2: Top: 5-run evaluation of test mean squared error on the fully-
synthetic (left) and semi-synthetic cases (case). CSI-3 was removed on the
right due to very high errors. Bottom: how errors change as training sizes are
increased, CSI-1 vs. GRU-2 (left: fully-synthetic, right: semi-synthetic).
Table 2: P-values from two-sample t-tests against the null hypothesis that
models perform no better on average than its counterpart. Significance at
$0.05$ indicated with one asterisk, $0.001$ with two.

| Full Synthetic | Semi-Synthetic Spotify
---|---|---
Model | CSI-VAE-2 | CSI-VAE-3 | GRU-0 | GRU-1 | GRU-2 | CSI-VAE-2 | CSI-VAE-3 | GRU-0 | GRU-1 | GRU-2
CSI-VAE-1 | $<0.001^{**}$ | $<0.001^{**}$ | $<0.001^{**}$ | $<0.001^{**}$ | $<0.001^{**}$ | $<0.001^{**}$ | $<0.001^{**}$ | $<0.001^{**}$ | $<0.001^{**}$ | $<0.001^{**}$
CSI-VAE-2 | $\backslash$ | $0.057$ | $<0.001^{**}$ | $<0.001^{**}$ | $0.249$ | $\backslash$ | $<0.001^{**}$ | $0.058$ | $0.081$ | $<0.05^{*}$
CSI-VAE-3 | $\backslash$ | $\backslash$ | $<0.001^{**}$ | $<0.001^{**}$ | $<0.05^{*}$ | $\backslash$ | $\backslash$ | $<0.001^{**}$ | $<0.001^{**}$ | $<0.001^{**}$
GRU-0 | $\backslash$ | $\backslash$ | $\backslash$ | $0.990$ | $<0.001^{**}$ | $\backslash$ | $\backslash$ | $\backslash$ | $0.705$ | $0.129$
GRU-1 | $\backslash$ | $\backslash$ | $\backslash$ | $\backslash$ | $<0.001^{**}$ | $\backslash$ | $\backslash$ | $\backslash$ | $\backslash$ | $0.061$

#### Datasets and oracular simulators.

The first step to assess intervention predictions is to build a set of proper
ground truth test beds, where we can control different levels of combinations
of interventions. We build two sets of oracular simulators. (1) The Fully-
synthetic simulator is constructed primarily based on random models following
our parameterization in Section 2. (2) The Semi-synthetic simulator is
constructed based on a real-world dataset from
Spotify333https://open.spotify.com/ which is aims to predict skip behaviour of
users based on its past interaction history. See details in Appendix C. For
each type of simulator, we generate $5$ different versions with different
random seeds. From simulator (1), we construct a simulated datasets of size
$50,000$, containing $5$ different interventions happening to the units at any
time after an initial burn-in period of $T_{0}=10$, although only maximum $3$
of these can be observed for any given unit, with $r=5$ being the true
dimensions of the model. For simulator (2), we construct simulated datasets of
size $3,000$, again containing $5$ different interventions, an initial period
$T_{0}=25$, a maximum 3 different interventions per unit and $r=10$. For
making predictions, we predict outcomes for interventions not applied yet
within any given unit (i.e. at least from the $2$ options left). In simulator
(2) parameters $\phi$ and $\beta$ are learned from real-world data.
Interventions are artificial, but inspired by the process of showing different
proportions of track types to an user in a Spotify-like environment. For both
setups, we use a data ratio of $0.7,0.1,0.2$ for training, validation and
test, respectively. We report the final results in Table 1 with another
constructed holdout set ($5,000$ for fully-synthetic and $1,000$ for semi-
synthetic).

#### Compared models.

We implemented three variations: (1) CSI-VAE-1 follows exactly our setup in
Section 2; (2) CSI-VAE-2 can be considered as an ablation study, which relaxes
the product form of Eq. 1 and replace it with a black-box MLP applied directly
to $(\phi_{n}^{t},\beta_{n},\psi_{n}^{t})$ (and hence may not guarantee
identifiability); (3) CSI-VAE-3 as another ablation study, where the equation
for $\psi$ (Eq. (2)) is replaced by a black-box function taking the sequence
of interventions $D$ as a standard time series. We compare our model against:
(1) GRU-0, a black-box gated recurrent unit (GRU, Cho et al. (2014)) composed
with a MLP, using only the past history of $X_{n}^{1:t}$ and $Z_{n}$; (2)
GRU-1, another GRU composed with a MLP that takes into account also the latest
intervention $D_{n}^{t}$; and (3) GRU-2, which uses not only $D_{n}^{t}$ but
also the entire past history $D_{n}^{1:t}$ just like CSI-VAE. In general,
GRU-2 can be considered as a very strong and generic black-box baseline model.

#### Results.

Each experiment was repeated 5 times, using Adam Kingma and Ba (2015) at a
learning rate of 0.01, with 50 epochs in the fully-synthetic case and 100 for
the semi-synthetic, which was enough for convergence. We selected the best
iteration point based on small holdout set. The main results are presented in
Tables 1 and 2, which shows the superiority of our model against the strong
baselines. We also observed that the identifiable results and compositional
interactions of intervention effect are both critical, as evidenced by the
drop in performance for CSI-VAE 2 and 3 in Table 1. In Appendix D, we provide
the following further experiments: (1) different choices of $r$ (to summarize:
using $r$ less than the true value gracefully underfits, while there is
evidence of small overfitting for choices of $r$ which overshoot the true
value but where we can enforce regularization techniques; (2) different size
of the training data (to summarize: even with more data provided, our model
consistently outperforms the blackbox models by a significant margin. We show
that a generic blackbox cannot solve this problem by simply feeding in more
data); and (3) a demonstration of conformal prediction that allow us to better
calibrate the predictive coverage compared to vanilla model-based prediction.

##  6 Conclusion

We introduced an approach for predicting potential outcomes with a careful
accounting of when extrapolation to unseen sequences are warranted. Lessons.
Embeddings are important given sparse categorical sequential data see e.g.
(Vaswani et al., 2017). However, large combinatorial interventional spaces
benefit from models that more carefully lay down conditions for
identification. Limitations. Unlike the traditional synthetic control
literature Abadie (2021), we assume a model for time effects based on
autoregression and truncated or parametric temporal-intervention interactions.
While it is possible to empirically evaluate the prediction abilities of the
model using a validation sample, high-stakes applications (such as major
interventions to counteract the effect of a pandemic) should take into
consideration that uncontrolled distribution shifts may take place, and
careful modeling of such shifts should be added to any analytical pipeline to
avoid damaging societal implications. Future Work. Besides allowing for an
explicit parameterization of drifts, an explicit account for modeling when
unmeasured confounding may have been presented in the selection of past
interventions is necessary to increase the applicability of the method to more
sources of data.

## Acknowledgements

We thank Zhenwen Dai and Georges Dupret for many helpful comments and
suggestions. JY and RS were partially funded by the EPSRC Open Fellowsship
EP/W024330/1. RS acknowledges further funding from the EPSRC AI Hub for
Causality in Healthcare AI with Real Data, EP/Y028856/1.

## References

  * Abadie [2021] A. Abadie.  Using synthetic controls: Feasibility, data requirements, and methodological aspects.  _Journal of Economic Literature_ , 59:391–425, 2021.
  * Agarwal et al. [2023a] A. Agarwal, A. Agarwal, and S. Vijaykumar.  Synthetic combinations: A causal inference framework for combinatorial interventions.  _Advances in Neural Information Processing Systems 36 (NeurIPS 2023)_ , pages 19195–19216, 2023a.
  * Agarwal et al. [2023b] A. Agarwal, M. Dahleh, D. Shah, and D. Shen.  Causal matrix completion.  _Proceedings of Thirty Sixth Conference on Learning Theory_ , pages 3821–3826, 2023b.
  * Agarwal et al. [2023c] A. Agarwal, D. Shah, and D. Shen.  Synthetic interventions.  _arXiv:2006.07691 (econ)_ , 2023c.
  * Aglietti et al. [2020] V. Aglietti, X. Lu, A. Paleyes, and J. Gonz’alez.  Causal Bayesian optimization.  _Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics_ , 34:3155–3164, 2020.
  * Angelopoulos et al. [2024] Anastasios Angelopoulos, Emmanuel Candes, and Ryan J Tibshirani.  Conformal pid control for time series prediction.  _Advances in Neural Information Processing Systems_ , 36, 2024.
  * Athey et al. [2021] S. Athey, M. Bayati, N. Doudchenko, G. Imbens, and K. Khosravi.  Matrix completion methods for causal panel data models.  _Journal of the American Statistical Association_ , 116:1716–1730, 2021.
  * Bravo-Hermsdorff et al. [2023] G. Bravo-Hermsdorff, D. Watson, J. Yu, J. Zeitler, and R. Silva.  Intervention generalization: A view from factor graph models.  _Advances in Neural Information Processing Systems 36 (NeurIPS 2023)_ , pages 43662–43675, 2023.
  * Brodersen et al. [2015] K. Brodersen, F. Gallusser, J. Koehler, N. Remy, and S. Scott.  Inferring causal impact using Bayesian structural time-series models.  _Annals of Applied Statistics, to appear_ , 2015.
  * Brost et al. [2019] Brian Brost, Rishabh Mehrotra, and Tristan Jehan.  The music streaming sessions dataset.  In _The World Wide Web Conference_ , pages 2594–2600, 2019.
  * Chakraborty and Moodie [2013] B. Chakraborty and E. Moodie.  _Statistical Methods for Dynamic Treatment Regimes: Reinforcement Learning, Causal Inference, and Personalized Medicine_.  Springer, 2013.
  * Chernozhukov et al. [2021] Victor Chernozhukov, Kaspar Wuthrich, and Yinchu Zhu.  An exact and robust conformal inference method for counterfactual and synthetic controls.  _Journal Of The American Statistical Association_ , 116(536):1849–1864, OCT 2 2021.  ISSN 0162-1459.  doi: 10.1080/01621459.2021.1920957.
  * Cho et al. [2014] K. Cho, B. van Merriënboer, C. Gulcehre, D. Bahdanau, F. Bougares, H. Schwenk, and Y. Bengio.  Learning phrase representations using RNN encoder–decoder for statistical machine translation.  _Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP 2014)_ , pages 1724–1734, 2014.
  * Chung et al. [2014] Junyoung Chung, Caglar Gulcehre, KyungHyun Cho, and Yoshua Bengio.  Empirical evaluation of gated recurrent neural networks on sequence modeling.  _arXiv preprint arXiv:1412.3555_ , 2014.
  * Correa and Bareinboim [2020] J. Correa and E. Bareinboim.  A calculus for stochastic interventions:causal effect identification and surrogate experiments.  _Proceedings of the AAAI Conference on Artificial Intelligence (AAAI 2020)_ , 34:10093–10100, 2020.
  * Dawid [2021] P. Dawid.  Decision-theoretic foundations for statistical causality.  _Journal of Causal Inference_ , 9(56):39–77, 2021.
  * Didelez et al. [2006] V. Didelez, A. P. Dawid, and S. Geneletti.  Direct and indirect effects of sequential treatments.  _Proceedings of the Twenty-Second Conference on Uncertainty in Artificial Intelligence (UAI2006)_ , 2006.
  * Evans and Didelez [2024] R. Evans and V. Didelez.  Parameterizing and simulating from causal models.  _Journal of the Royal Statistical Society. Series B (Statistical Methodology)_ , 2024.  URL https://doi.org/10.1093/jrsssb/qkad058.
  * Fontana et al. [2023] Matteo Fontana, Gianlca Zeni, and Simone Vantini.  Conformal prediction: A unified review of theory and new challenges.  _Bernoulli_ , 29(1):1–23, FEB 2023.  ISSN 1350-7265.  doi: 10.3150/21-BEJ1447.
  * Gibbs and Candes [2021] Isaac Gibbs and Emmanuel Candes.  Adaptive conformal inference under distribution shift.  _Advances in Neural Information Processing Systems_ , 34:1660–1672, 2021.
  * Gui et al. [2023] Y. Gui, R. Barber, and C. Ma.  Conformalized matrix completion.  _Advances in Neural Information Processing Systems_ , 36:4820–4844, 2023.
  * Hernán and J [2020] M. Hernán and Robins J.  _Causal Inference: What If_.  Chapman $\&$ Hall/CRC, 2020.
  * Hochreiter and Schmidhuber [1997] S. Hochreiter and J. Schmidhuber.  Long short-term memory.  _Neural Computation_ , 9:1735–1780, 1997.
  * Imbens and Rubin [2015] G. W. Imbens and D. B. Rubin.  _Causal Inference in Statistics, Social, and Biomedical Sciences: An Introduction_.  Cambridge University Press, 2015.
  * Kaddour et al. [2021] J. Kaddour, Y. Zhu, Q. Liu, M. J. Kusner, and R. Silva.  Causal effect inference for structured treatments.  _Advances in Neural Information Processing Systems 34 (NeurIPS 2021)_ , pages 24841–24854, 2021.
  * Kingma and Ba [2015] D. Kingma and J. Ba.  Adam: A method for stochastic optimization.  _3rd International Conference for Learning Representations,_ , 2015.
  * Kingma and Welling [2013] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  _arXiv preprint arXiv:1312.6114_ , 2013.
  * Lei and Candès [2021] L. Lei and E. J. Candès.  Conformal inference of counterfactuals and individual treatment effects.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 83(5):911–938, 2021.
  * Mnih and Gregor [2014] Andriy Mnih and Karol Gregor.  Neural variational inference and learning in belief networks.  In _International Conference on Machine Learning_ , pages 1791–1799. PMLR, 2014.
  * Murphy [2003] S. Murphy.  Optimal dynamic treatment regimes.  _Journal of the Royal Statistical Society Series B_ , 6(2):331–355, 2003.
  * Nabi et al. [2022] R. Nabi, T. McNutt, and I. Shpitser.  Semiparametric causal sufficient dimension reduction of multidimensional treatments.  _The 38th Conference on Uncertainty in Artificial Intelligence_ , pages 1445–1455, 2022.
  * Pearl [2009] J. Pearl.  _Causality: Models, Reasoning and Inference, 2nd edition_.  Cambridge University Press, 2009.
  * Rezende et al. [2014] Danilo Jimenez Rezende, Shakir Mohamed, and Daan Wierstra.  Stochastic backpropagation and approximate inference in deep generative models.  In _International conference on machine learning_ , pages 1278–1286. PMLR, 2014.
  * Robins [1986] James Robins.  A new approach to causal inference in mortality studies with a sustained exposure period—application to control of the healthy worker survivor effect.  _Mathematical modelling_ , 7(9-12):1393–1512, 1986.
  * Robins [1997] James M Robins.  Causal inference from complex longitudinal data.  In _Latent variable modeling and applications to causality_ , pages 69–117. Springer, 1997.
  * Saengkyongam et al. [2024] S. Saengkyongam, Rosenfeld, P. Ravikumar, N. Pfister, and J. Peters.  Identifying representations for intervention extrapolation.  _International Conference on Learning Representations (ICLR 2024)_ , 2024.
  * Saito et al. [2023] Y. Saito, R. Qingyang, and T. Joachims.  Off-policy evaluation for large action spaces via conjunct effect modeling.  _Proceedings of the 39th International Conference on Machine Learning (ICML 2023)_ , pages 29734––29759, 2023.
  * Saito et al. [2024] Y. Saito, J. Yao, and T. Joachims.  POTEC: Off-policy learning for large action spaces via two-stage policy decomposition.  _Proceeding of the 40th International Conference on Machine Learning (ICML 2024)_ , 2024.
  * Schäfer and Zimmermann [2007] A.M. Schäfer and H.-G. Zimmermann.  Recurrent neural networks are universal approximators.  _International Journal of Neural Systems_ , 17(04):253–263, 2007.
  * Spirtes et al. [2000] P. Spirtes, C. Glymour, and R. Scheines.  _Causation, Prediction and Search_.  Cambridge University Press, 2000.
  * Squires et al. [2022] C. Squires, D. Shen, A. Agarwal, D. Shah, and C. Uhler.  Causal imputation via synthetic interventions.  _First Conference on Causal Learning and Reasoning (CLeaR 2021)_ , pages 688–711, 2022.
  * Sutton and Barto [2018] R. Sutton and A. Barto.  _Reinforcement Learning: An Introduction_.  MIT Prsss, 2018.
  * Tibshirani et al. [2019] Ryan J Tibshirani, Rina Foygel Barber, Emmanuel Candes, and Aaditya Ramdas.  Conformal prediction under covariate shift.  _Advances in neural information processing systems_ , 32, 2019.
  * Vaswani et al. [2017] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. Gomez, Ł. Kaiser, and I. Polosukhin.  Attention is all you need.  _Advances in Neural Information Processing Systems 31 (NIPS 2017)_ , pages 6000–6010, 2017.
  * Vovk et al. [2005] V. Vovk, A. Gammerman, and G. Shafer.  _Algorithmic Learning in a Random World_.  Springer, 2005.
  * Yao et al. [2019] L. Yao, S. Li, Y. Li, H. Xue, J. Gao, and A. Zhang.  On the estimation of treatment effect with text covariates.  _Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence (IJCAI-19)_ , pages 4106–4113, 2019.
  * Zeng et al. [2024] Z. Zeng, S. Balakrishnan, Y. Han, and E. Kennedy.  Causal inference with high-dimensional discrete covariates.  _arXiv:2405.00118 (math)_ , 2024.
  * Zhang and Bareinboim [2020] J. Zhang and E. Bareinboim.  Designing optimal dynamic treatment regimes: A causal reinforcement learning approach.  _Proceedings of the 37th International Conference on Machine Learning (ICML 2020)_ , 2020.
  * Zhang et al. [2023] J. Zhang, K. Greenewald, C. Squires, A. Srivastava, K. Shanmugam, and C. Uhler.  Identifiability guarantees for causal disentanglement from soft interventions.  _Advances in Neural Information Processing Systems_ , 36:50254–50292, 2023.

## Appendix

##  Appendix A General Background

Much of traditional causal inference aimed at estimating an outcome of a fixed
treatment that does not vary over time Imbens and Rubin [2015]. However, many
causal questions involves treatments that vary over time. For example, we
might be interested in estimating the causal effects of medical treatments,
lifestyle habits, employment status, marital status, occupation exposures,
etc. Because these treatments may take different values for a single
individual over time, hence they are often refereed to as time-varying
treatments [Hernán and J, 2020].

In many designs, it is assumed that a time-fixed treatment applied to
individuals in our population happens at the same time (i.e. there is a shared
global clock for everyone). A more relaxed setting (and more practical, as
described in this paper) is that (categorical) interventions can happen to
different individuals at different time steps, potentially resulting in no two
identically distributed time-series. The purpose now is to estimate the causal
effect or predict expected outcomes under time-varying treatments as a time-
series, rather than estimating missing potential outcome snapshots at fixed
time interval that has passed already.

One typical way of designing experiments under a binary treatment ($0$ or $1$)
[Robins, 1986, 1997] of time-varying treatments is to have
$\bar{d}_{0}=(0,0,\cdots,0)$ (never treated) and $\bar{d}_{1}=(1,1,\cdots,1)$
(always treated) as the two options in the design space. This allows us to
answer what the average expected outcome contrasted between these two options.
A fine-grained estimation at every time step will require further levels of
contrast under two regimes, such as $d=(0,1,\cdots,1)$ and $d=(1,0,\cdots,0)$.
However, the number of combinations grows exponentially on $T$, the number of
time steps, even when only binary treatments are at the stake. In this paper,
we focus on static treatment strategies where the treatment is randomly
assigned (also referred to as a sequential randomized experiment), although
the next intervention in a sequential randomised experiment can depend on
previous step intervention, as long as there is no influence from unobserved
confounders.

Conditional ignorability holds when we are able to block unmeasured
confounders using observable confounders only. Similarly, in a time-varying
setting, we achieve sequential conditional ignorability by conditioning on the
time-varying covariates, including past actions (Hernán and J [2020], and
Chapter 4 of Pearl [2009]). In our manuscript, we will have little to say
about the uses of conditional sequential ignorability and other adjustment
techniques for observational studies, as these ideas are well-understood and
can be applied on top of our main results. Hence, we decided to focus on
treatments that act as-if randomized so that we can explore in more detail our
contributions.

##  Appendix B Proofs

#### Proposition 1

_Let
$f(n,d^{1:T},x^{1:T},z):=f_{d}^{\mathsf{T}}(d^{1:T})f_{nxz}(n,x^{1:T},z)$,
where function sequences $f_{d}(d^{1:T})$ and $f_{nxz}(n,x^{1:T},z)$ are
defined for all $T\in\mathbb{N}^{+}$ and have codomain $\mathbb{R}^{r_{g}}$.
Furthermore, define $f_{dl}(d^{1:T})$ as in Eq. (4) with $r_{h}=1$, and the
$l$-th entry of $f_{nxz}(n,x^{1:T},z)$ as
$f_{nxzl}(n,x^{1:T},z):=u_{l}(n)^{\mathsf{T}}v_{l}(x^{1:T},z)$, each of the
latter having codomain $\mathbb{R}^{r_{s^{\prime}}}$ and where
$v_{l}(\cdot,\cdot)$ is defined recursively for any $T$, similarly
$g_{l}(\cdot,\cdot,\cdot)$. Then there exists some integer $r$ and three
functions $a,b$ and $c$ with codomain $\mathbb{R}^{r}$ so that
$f(n,d^{1:T},x^{1:T},z)=\sum_{l=1}^{r}a_{l}(n)\times b_{l}(x^{1:T},z)\times
c_{l}(d^{1:T},T)$. Moreover $c_{l}(d^{1:T},T)=\prod_{t=1}^{T}m_{l}(d^{t},t,T)$
for functions $m_{l}:\mathbb{N}^{3}\rightarrow\mathbb{R}$._

#### Proof.

The proof is tedious and the result intuitive, but they help to formalize the
main motivation for Eqs. (1) and (3).

$\begin{array}[]{rcl}f_{d}^{\mathsf{T}}(d^{1:T})f_{nxz}(n,x^{1:T},z)&=&\sum_{l^{\prime}=1}^{r_{s}}f_{dl^{\prime}}(d^{1:T})f_{nxzl^{\prime}}(n,x^{1:T},z)\\\
\\\
&=&\sum_{l^{\prime}=1}^{r_{s}}f_{dl^{\prime}}(d^{1:T})\sum_{l^{\prime\prime}}^{r_{s^{\prime}}}u_{l^{\prime}l^{\prime\prime}}(n)v_{l^{\prime}l^{\prime\prime}}(x^{1:T},z)\\\
\\\
&=&\sum_{l^{\prime\prime}=1}^{r_{s^{\prime}}}\sum_{l^{\prime}=1}^{r_{s}}u_{l^{\prime}l^{\prime\prime}}(n)v_{l^{\prime}l^{\prime\prime}}(x^{1:T},z)f_{dl^{\prime}}(d^{1:T})\\\
\\\ &=&\sum_{l=1}^{r}a_{l}(n)b_{l}(x^{1:T},z)c_{l}(d^{1:T}),\end{array}$

where $r=r_{s}\times r_{s^{\prime}}$ and
$a_{l}(n):=u_{l^{\prime}l^{\prime\prime}}(n)$ for $l=l^{\prime}\times
l^{\prime\prime}$;
$b_{l}(x^{1:T},z):=v_{l^{\prime}l^{\prime\prime}}(x^{1:T},z)$ for
$l=l^{\prime}\times l^{\prime\prime}$; and
$c_{l}(d^{1:T})=f_{dl^{\prime}}(d^{1:T})$ for $l^{\prime}=\lceil
l/r_{s^{\prime}}\rceil$, where $\lceil\cdot\rceil$ is the ceiling function.

Finally, each $c_{l}(d^{1:T})$ is by definition some
$f_{dl^{\prime}}(d^{1:T})$. From Eq. (4) and the assumptions,
$f_{dl^{\prime}}(d^{1:T})=g_{l^{\prime}}(d^{1:T},T,T)$, and
$g_{l^{\prime}}(d^{1:t},t^{\prime},t):=h_{l^{\prime}}(g_{l}(d^{1:t^{\prime}-1},t^{\prime}-1,t),d^{t^{\prime}},t^{\prime},t)$
for $t^{\prime}>0$ and $d>0$, with otherwise
$h_{l^{\prime}_{1}}(g_{l^{\prime}}(d^{1:t},t^{\prime}-1,t))h_{l^{\prime}_{2}}(d^{t},t^{\prime},t)$.
For $t^{\prime}=0$, $h_{l^{\prime}}(v,d,t^{\prime},t)\equiv 1$, and for
$t^{\prime}>0$ and $d=0$, $h_{l^{\prime}}(v,d,t^{\prime},t)\equiv
g_{l^{\prime}}(d^{1:t^{\prime}-1},t^{\prime}-1,t)$. By recursive application,
this implies
$f_{dl^{\prime}}(d^{1:T})=\prod_{t^{\prime}=1:d^{t^{\prime}}>0}^{T}h_{l^{\prime}_{2}}(d^{t^{\prime}},t^{\prime},T)$
and the result follows. $\Box$

#### Proposition 2

_Suppose we have a measurable function
$f:\bigcup_{t=2}^{T}\mathcal{X}^{1:t-1}\times\bigcup_{t=1}^{T}\mathcal{D}^{t}\times\mathcal{Z}\rightarrow\mathbb{R}^{d}$,
$t\leq T$, where $\mathcal{X}^{1:t-1}$ and $\mathcal{Z}$ are compact,
$\mathcal{D}^{t}$ is finite. For $d=1,\cdots,d_{max}$, let
$m_{d}:\bigcup_{t=1}^{T}\mathcal{D}^{t}\rightarrow\\{1,\cdots,t\\}$ with
$m_{d}(\mathbf{d})=t\text{ if }\exists t\;s.t.\;\mathbf{d}_{t}=d\;\text{else
}0$. There exist $L=O(T^{d_{max}+1})$ such that for any $\epsilon$, there are
measurable functions
$\phi_{l}:\mathcal{X}^{1:t-1}\times\mathcal{Z}\rightarrow\mathbb{R}^{d}$ and
$\psi_{ld^{\prime}}:\mathcal{D}^{t}\rightarrow\mathbb{R}$ and real numbers
$\beta_{l}$ for $l=1,\cdots,r$ and $d=1,\cdots,d_{max}$,
$\hat{f}(\mathbf{x},\mathbf{d},z)=\sum_{l=1}^{L}\phi_{l}(\mathbf{x},z)\beta_{l}\prod_{d=1}^{d_{max}}\psi_{ld^{\prime}}(m_{d^{\prime}}(\mathbf{d}))$
such that_

$\displaystyle\text{sup}_{\mathbf{x}\in\bigcup_{t=2}^{T}\mathcal{X}^{1:t-1},\mathbf{d}\in\bigcup_{t=1}^{T}\mathcal{D}^{t},z\in\mathcal{Z}}\left|f(\mathbf{x},\mathbf{d},z)-\hat{f}(\mathbf{x},\mathbf{d},z)\right|\leq\epsilon$
(11)

_Proof._ First note that we can construct an equivalent definition of $D^{t}$
for some fixed $t$ as follows. We for now call the new object $\bar{D}^{t}$.

$\bar{\mathcal{D}}^{t}=\left\\{\mathbf{v}\in\mathbb{R}^{d_{max}}\bigg{|}v_{d}\in\\{0,\cdots,t\\}\forall
d;\;v_{d^{\prime}}\neq v_{d}\text{ whenever $d^{\prime}\neq d$ and $v_{d}\neq
0$}\right\\}$ (12)

Note that there is a one-to-one correspondence between $\bar{\mathcal{D}}^{t}$
and $\mathcal{D}^{t}$, and denote the bijection connecting them by
$g:\bar{\mathcal{D}}^{t}\rightarrow\mathcal{D}^{t}$. Since both sets are
finite, we can approximate $f$ with $\hat{f}$, if and only if we can
approximate $f\circ i\times g\times i$ by $\hat{f}\circ i\times g\times i$. So
wlog just consider $\bar{\mathcal{D}}^{t}$ as $\mathcal{D}^{t}$.

The set of all sequences of treatments is the union of $\mathcal{D}^{t}$ over
$t=1,...,T$, for some known $T$:

$\displaystyle\mathcal{D}^{\mathcal{T}}$
$\displaystyle=\bigcup_{t=1}^{T}\mathcal{D}^{t}\;\;\text{known $T$.}$ (13)

Note that

$\displaystyle|\mathcal{D}^{t}|$ $\displaystyle<(t+1)^{d_{max}}$ (14)
$\displaystyle|\mathcal{D^{T}}|$
$\displaystyle<\sum_{t=1}^{T}(t+1)^{d_{max}}\;\;$ (15)

therefore the size of $\mathcal{D^{T}}$ is polynomial in $T$.

Now note, fixing the ground truth function $f$ at $D\in\mathcal{D^{T}}$ gives
me a new function
$\bigcup_{t=1}^{T}\mathcal{X^{1:t-1}}\times\mathcal{Z}\rightarrow\mathcal{X}$.
Analogously, fixing $\hat{f}$ at some $D\in\mathcal{D^{T}}$ also give a new
function in the same function space. Since $\mathcal{D^{T}}$ is finite, we can
enumerate its elements. Thus for every $D^{(i)}\in\mathcal{D^{T}}$,
$i=1,\cdots,|\mathcal{D^{T}}|$, write down $\hat{f}$ fixed at $D^{(i)}$:

$\displaystyle\hat{f}(X^{1:t-1},D^{(i)},Z)$
$\displaystyle=\sum_{l=1}^{r}\phi_{l}(X^{1:t-1},Z)\beta_{l}\psi_{l1}(D^{(i)}_{1})\psi_{l2}(D^{(i)}_{2})\cdots\psi_{ld_{max}}(D^{(i)}_{d_{max}})$
(16)

Then we can express the vector field:

$\displaystyle\begin{pmatrix}\hat{f}(\cdot,D^{(1)},\cdot)\\\ \vdots\\\
\hat{f}(\cdot,D^{(|\mathcal{D^{T}}|)},\cdot)\end{pmatrix}$
$\displaystyle=\underbrace{\begin{pmatrix}\beta_{1}\prod_{d=1}^{d_{max}}\psi_{11}(D^{(1)}_{1})&\cdots&\beta_{r}\prod_{d=1}^{d_{max}}\psi_{r1}(D^{(1)}_{1})\\\
\vdots&\vdots&\vdots\\\
\beta_{1}\prod_{d=1}^{d_{max}}\psi_{11}(D^{(|\mathcal{D^{T}}|)}_{1})&\cdots&\beta_{r}\prod_{d=1}^{d_{max}}\psi_{r1}(D^{(|\mathcal{D^{T}}|)}_{1})\end{pmatrix}}_{W}\begin{pmatrix}\phi_{1}(\cdot,\cdot)\\\
\vdots\\\ \phi_{r}(\cdot,\cdot)\end{pmatrix}$ (17)

By choosing $r=|\mathcal{D^{T}}|$ and appropriate values for $\beta_{l}$ and
$\psi_{ld}(D^{(i)})$, we can make $W$ invertible.

Now consider the vector field given by the ground truth function fixed at all
values of $D\in\mathcal{D^{T}}$:

$\displaystyle\begin{pmatrix}f(\cdot,D^{(1)},\cdot)\\\ \vdots\\\
f(\cdot,D^{(|\mathcal{D^{T}}|)},\cdot)\end{pmatrix}$ (18)

And consider its transformation under $W^{-1}$:

$\displaystyle W^{-1}\begin{pmatrix}f(\cdot,D^{(1)},\cdot)\\\ \vdots\\\
f(\cdot,D^{(|\mathcal{D^{T}}|)},\cdot)\end{pmatrix}$ (19)

Since recurrent neural networks are universal approximators Schäfer and
Zimmermann [2007], for every $\delta>0$ we can choose
$\phi_{1},\cdots,\phi_{r}$ so that

$\displaystyle\text{sup}_{\mathbf{x}\in\bigcup_{t=2}^{T}\mathcal{X}^{1:t-1},\;z\in\mathcal{Z}}\left\|\begin{pmatrix}\phi_{1}(\mathbf{x},z)\\\
\vdots\\\
\phi_{r}(\mathbf{x},z)\end{pmatrix}-W^{-1}\begin{pmatrix}f(\mathbf{x},D^{(1)},z)\\\
\vdots\\\
f(\mathbf{x},D^{(|\mathcal{D^{T}}|)},z)\end{pmatrix}\right\|^{2}_{2}\leq\delta$
(20)

Since $W$ is a finite matrix, its operator norm is well-defined. Then

$\displaystyle\text{sup}_{\mathbf{x}\in\bigcup_{t=2}^{T}\mathcal{X}^{1:t-1},\;z\in\mathcal{Z}}\left\|W\begin{pmatrix}\phi_{1}(\mathbf{x},z)\\\
\vdots\\\
\phi_{r}(\mathbf{x},z)\end{pmatrix}-\begin{pmatrix}f(\mathbf{x},D^{(1)},z)\\\
\vdots\\\
f(\mathbf{x},D^{(|\mathcal{D^{T}}|)},z)\end{pmatrix}\right\|^{2}_{2}$ (21)
$\displaystyle\leq$
$\displaystyle\text{sup}_{\mathbf{x}\in\bigcup_{t=2}^{T}\mathcal{X}^{1:t-1},\;z\in\mathcal{Z}}\|W\|^{2}_{op}\left\|\begin{pmatrix}\phi_{1}(\mathbf{x},z)\\\
\vdots\\\
\phi_{r}(\mathbf{x},z)\end{pmatrix}-W^{-1}\begin{pmatrix}f(\mathbf{x},D^{(1)},z)\\\
\vdots\\\
f(\mathbf{x},D^{(|\mathcal{D^{T}}|)},z)\end{pmatrix}\right\|^{2}_{2}$ (22)
$\displaystyle\leq$
$\displaystyle\|W\|^{2}_{op}\text{sup}_{\mathbf{x}\in\bigcup_{t=2}^{T}\mathcal{X}^{1:t-1},\;z\in\mathcal{Z}}\left\|\begin{pmatrix}\phi_{1}(\mathbf{x},z)\\\
\vdots\\\
\phi_{r}(\mathbf{x},z)\end{pmatrix}-W^{-1}\begin{pmatrix}f(\mathbf{x},D^{(1)},z)\\\
\vdots\\\ f(\mathbf{x},D^{(|\mathcal{D^{T}}|)},z)\end{pmatrix}\right\|_{2}$
(23) $\displaystyle\leq$ $\displaystyle\|W\|^{2}_{op}\delta$ (24)

But

$\displaystyle\text{sup}_{\mathbf{x}\in\bigcup_{t=2}^{T}\mathcal{X}^{1:t-1},\;\mathbf{d}\in\mathcal{D^{T}}\;z\in\mathcal{Z}}|f(\mathbf{x},\mathbf{d},z)-\hat{f}(\mathbf{x},\mathbf{d},z)|$
(25) $\displaystyle=$
$\displaystyle\left(\text{sup}_{\mathbf{x}\in\bigcup_{t=2}^{T}\mathcal{X}^{1:t-1},\;\;z\in\mathcal{Z}}\text{max}_{\mathbf{d}\in\mathcal{D^{T}}}|f(\mathbf{x},\mathbf{d},z)-\hat{f}(\mathbf{x},\mathbf{d},z)|^{2}\right)^{1/2}$
(26) $\displaystyle\leq$
$\displaystyle\left(\text{sup}_{\mathbf{x}\in\bigcup_{t=2}^{T}\mathcal{X}^{1:t-1},\;\;z\in\mathcal{Z}}\sum_{\mathbf{d}\in\mathcal{D^{T}}}|f(\mathbf{x},\mathbf{d},z)-\hat{f}(\mathbf{x},\mathbf{d},z)|^{2}\right)^{1/2}$
(27) $\displaystyle=$
$\displaystyle\left(\text{sup}_{\mathbf{x}\in\bigcup_{t=2}^{T}\mathcal{X}^{1:t-1},\;\;z\in\mathcal{Z}}\left\|W\begin{pmatrix}\phi_{1}(\mathbf{x},z)\\\
\vdots\\\
\phi_{r}(\mathbf{x},z)\end{pmatrix}-\begin{pmatrix}f(\mathbf{x},D^{(1)},z)\\\
\vdots\\\
f(\mathbf{x},D^{(|\mathcal{D^{T}}|)},z)\end{pmatrix}\right\|^{2}_{2}\right)^{1/2}$
(28) $\displaystyle\leq$ $\displaystyle\|W\|^{2}_{op}\delta$ (29)

So choosing $\delta=\epsilon/\|W\|^{2}_{op}$ would give us the desired result.

$\Box$

#### Proposition 3

_Let $b_{n}$ be a $(T_{0_{n}}-1)\times 1$ vector, where each entry $t$ is
given by
$\mathbb{E}[X_{n}^{t+1}(d_{n}^{1:{t+1}})~{}|~{}x_{n}^{1:t}(d_{n}^{1:t}),z_{n}]$.
Under Assumptions 1 and 2, $\beta_{n}$ is identifiable from $A_{n}$ and
$b_{n}$._

#### Proof.

Let $b_{n}$ be a $(T_{0_{n}}-1)\times 1$ vector where each entry $t$ is given
by $\mathbb{E}[X_{n}^{t+1}(0^{1:t+1})~{}|~{}x_{n}^{1:t}(0^{1:t}),z_{n}]$. We
know from basic least-squares theory, and the given assumptions, that as long
as the number $T_{0}\geq r$ of observations is larger than or equal to the
dimension $r$, the system is (over-)determined and $\beta_{n}=A_{n}^{+}b$.
$\Box$

#### Theorem 1

_Assume we have a dataset of $N$ units and $T$ time points
$(d_{n}^{1:T},x_{n}^{1:T},z_{n})$ generated by a model partially specified by
Eq. (1). Assume also knowledge of the conditional expectations
$\mathbb{E}[X_{n}^{t}(d_{n}^{1:t})~{}|~{}x_{n}^{t-1}(d_{n}^{1:t-1}),z_{n}]$
and basis functions $\phi_{l}(\cdot,\cdot)$ for all $l=1,2,\dots,r$ and all
$1\leq n\leq N$ and $1\leq t\leq T$. Then, under Assumptions 1-4 applied to
all individuals and all intervention levels $d$ that appear in our dataset, we
have that all $\beta_{n}$ and all $\psi_{l}(d,\cdot,\cdot)$ are identifiable._

#### Proof.

By Assumptions 1 and 2, and Proposition 3, all $\beta_{n}$ are identifiable.

For the next step, let

$\eta_{nl}^{t}:=\beta_{nl}\times\prod_{t^{\prime}=1}^{t}\psi_{l}(d_{n}^{t^{\prime}},t^{\prime},t),$

and

$\alpha_{dl}^{\Delta_{t}}:=\prod_{t^{\prime}=1}^{\Delta_{t}}\psi_{l}(d,1,t^{\prime}).$

Consider now the first time $t$ an intervention level $d>0$ is assigned to
someone in the entire dataset. By Assumption 3, there exists a dataset
$\mathcal{N}_{d}$ with the given properties where, for all
$n\in\mathcal{N}_{d}$, we have
$\eta_{nl}^{t}=\beta_{nl}\times\prod_{t^{\prime}=1}^{t}\psi_{l}(d_{n}^{t^{\prime}},t^{\prime},t)=\beta_{nl}$,
since all interventions prior to $t$ have been at level 0.

Let $b_{dt^{\prime}}$ be a $N_{d}\times 1$ vector where each entry $i$ is
given by
$\mathbb{E}[X_{n_{i}}^{t+t^{\prime}}(d_{n_{i}}^{1:t+t^{\prime}})~{}|~{}x_{n_{i}}^{1:t+t^{\prime}}(d_{n_{i}}^{1:t+t^{\prime}}),z_{n_{i}}]$,
for $t^{\prime}=1,2,\dots,k_{d}-1$. Consider the time-bounded case with free
parameters $w_{dl,1}$, $l=1,2,\dots,r$. Given Assumption 4, we can solve for
$\eta_{n}^{t+1}=A_{d_{1}}^{+}b_{d_{1}}$, the vector with entries
$\eta_{nl}^{t+1}$.

As
$\eta_{nl}^{t+1}=\eta_{nl}^{t}\times\psi_{l}(d,t+1,t+1)=\eta_{nl}^{t+1}\times
w_{dl,1}$, with non-zero $\eta_{nl}$ by Assumption 3, this also identifies
$w_{dl,1}$. Also by Assumption 3, at time point $t+2$ we have that
$\eta_{nl}^{t+1}=\eta_{nl}^{t}\times\psi_{l}(d,t+1,t+2)$ for all units in
$\mathcal{N}_{d}$, which by analogy to the previous paragraph, identifies
$\psi_{l}(d,t+1,t+2)$ and consequently $w_{dl,2}$. As this goes all the way to
$t+k_{d}-1$, this identifies all of $\psi_{l}(d,\cdots,\cdot)$.

If intervention $d$ level is the time-unbounded model of Eq. (3), by
Assumption 3, we also identify $\alpha_{d}^{1}$, $\alpha_{d}^{2}$ and
$\alpha_{d}^{3}$, as they are defined by functions
$\psi_{l}(d,t+1,t+t^{\prime})=\psi_{l}(d,1,t^{\prime})$ for $t^{\prime}=1,2,3$
that we established as identified by the reasoning in the previous paragraph.
As $\sigma(w_{1dl})\times w_{2dl}-\sigma(w_{1dl})^{3}\times
w_{2dl}=\alpha_{dl}^{1}-\alpha_{dl}^{3}$ and

$\displaystyle\frac{\sigma(w_{1dl})\times w_{2dl}-\sigma(w_{1dl})^{2}\times
w_{2dl}}{\sigma(w_{1dl})^{2}\times w_{2dl}-\sigma(w_{1dl})^{3}\times
w_{2dl}}=\frac{\alpha_{dl}^{1}-\alpha_{dl}^{2}}{\alpha_{dl}^{2}-\alpha_{dl}^{3}},$

solving for the above results in

$\begin{array}[]{rcl}w_{1dl}&=&\displaystyle\sigma^{-1}\left(\frac{\alpha_{dl}^{2}-\alpha_{dl}^{3}}{\alpha_{dl}^{1}-\alpha_{dl}^{2}}\right),\\\
\\\
w_{2dl}&=&\displaystyle\frac{\alpha_{dl}^{1}-\alpha_{dl}^{3}}{\sigma(w_{1dl})-\sigma(w_{1dl})^{3}},\\\
\\\ w_{3dl}&=&\alpha_{dl}^{1}-\sigma(w_{1dl})\times w_{2dl}.\end{array}$

So far we have established identification of the $\psi(d,\cdot,\cdot)$
functions for the first intervention level $d$ that appears in the dataset
(the value $d$ is not unique, as it is possible to assign different
intervention levels $d^{\prime}$ in parallel at time $t$, but to different
units). We now must show that, for the remaining interventions $d>0$ levels
take place in our dataset after time $t$.

Assume the induction hypothesis that $t_{s}$ is the $s$-th unique time an
intervention level $d>0$ is assigned to some unit in the dataset, and all
intervention levels $d^{\star}$ that appeared up prior to $t_{d}$ have their
functions $\psi_{l}(d^{\star},\cdot,\cdot)$ identified. We showed this is the
case for $t_{1}=t$, which is our base case. We assume this to be true for
$t_{s}$ and we will show that this also holds for $t_{s+1}$.

Let $d$ be an intervention level assigned at $t_{s+1}$ and let
$\mathcal{N}_{d}$ the corresponding subset of data assumed to exist by
Assumption 3. For any $n\in\mathcal{N}_{d}$, let $s_{-}$ be the last time a
non-zero intervention $d_{-}$ has been assigned to $n$ prior to $s+1$. Also by
Assumption 3, this must have taken place at least $k_{d_{-}}$ steps in the
past. By the induction hypothesis, $\psi_{l}(k_{d_{-}},\cdot,\cdot)$ is
identifiable, and therefore so is the case for all
$\psi_{l}(k_{d_{\prec}},\cdot,\cdot)$ for intervention levels $d_{\prec}$
taking place for the first time prior to $t_{s}$. This also implies that
$\eta_{nl}^{t_{s}}$ is identified, and the rest of the argument proceeds as in
the base case. $\Box$

#### Theorem 2

_Assume we have $N$ calibration samples,_

$\displaystyle
S_{n}^{t^{\prime}}=|X^{t^{\prime}}_{n}-f(D^{1:t_{n}}_{n},X^{1:t_{n}}_{n},Z_{n})|,\quad
t^{\prime}=t_{n}+\Delta<t,\quad n=1,\dots,N$ (30)

_where $t_{n}$ is the time user $n$ experienced the last intervention before
$t$. Assume there exists $\epsilon>0$ such that, for all $n$,_

$\displaystyle
p_{T}(S_{n}^{t+\Delta})=(1-\epsilon)p_{C}(S_{n}^{t+\Delta})+\epsilon
p_{\delta}(S_{n}^{t+\Delta}),$ (31)

_where $p_{T}$ and $p_{C}$ are the (unknown) densities of the test and
calibration distributions, $p_{\delta}$ is a bounded arbitrary shift density,
and $p_{min}=\min_{n=1,\dots,N}p_{C}(S_{n}^{t^{\prime}})>0$. Then,_

$\displaystyle{\rm Prob}\left(X_{N+1}^{t+\Delta}\in C\right)\geq
1-\alpha-\frac{1}{p_{min}}\frac{\epsilon}{1-\epsilon}.$ (32)

#### Proof.

Let $P_{T}$ and $P_{C}$ be the joint distributions of
$(X_{n}^{t^{\prime}},D_{n}^{t^{\prime}})$ for $t^{\prime}>t$ and
$t^{\prime}\leq t$ and
$w_{n}^{t^{\prime}}=\frac{p_{T}(X^{t^{\prime}}_{n},D^{t^{\prime}}_{n})}{p_{C}(X^{t^{\prime}}_{n},D^{t^{\prime}}_{n})}$.
For any $t^{\prime}>t$, we can estimate $P_{T}$ by weighting the calibration
samples with $w_{n}^{t^{\prime}}$, $n=1,\dots,N$, i.e.

$\displaystyle P_{T}(X,D)\approx Z^{-1}\sum_{n=1}^{N}w_{n}^{t^{\prime}}{\bf
1}((X,D)=(X_{n}^{t^{\prime}},D_{n}^{t^{\prime}})),\quad
Z=\sum_{n=1}^{N}w_{n}^{t^{\prime}}$ (33)

Conditional on the calibration samples, the $(1-\alpha)$-th empirical quantile
of $P_{T}(X,D)$ is

$\displaystyle Q_{T}(1-\alpha)={\rm
inf}_{q}\left\\{\left(Z^{-1}\sum_{n=1}^{N}w_{n}^{t^{\prime}}{\bf 1}(S_{n}\leq
q)\right)\geq 1-\alpha\right\\}$ (34) $\displaystyle={\rm
inf}_{q}\left\\{\left(\sum_{n=1}^{N}\left(\frac{w_{n}^{t^{\prime}}}{Z}-\frac{1}{N}+\frac{1}{N}\right){\bf
1}(S_{n}\leq q)\right)\geq 1-\alpha\right\\}$ (35) $\displaystyle={\rm
inf}_{q}\left\\{\left(\frac{1}{N}\sum_{n=1}^{N}{\bf 1}(S_{n}\leq q)\right)\geq
1-\alpha-\sum_{n=1}^{N}\left(\frac{w_{n}^{t^{\prime}}}{Z}-\frac{1}{N}\right){\bf
1}(S_{n}\leq q)\right\\}$ (36) $\displaystyle\geq{\rm
inf}_{q}\left\\{\left(\frac{1}{N}\sum_{n=1}^{N}{\bf 1}(S_{n}\leq q)\right)\geq
1-\alpha-\sqrt{N}\sqrt{\sum_{n=1}^{N}\left(\frac{w_{n}^{t^{\prime}}}{Z}-\frac{1}{N}\right)^{2}}\right\\}$
(37) $\displaystyle=Q_{C}(1-\alpha-{\rm gap})$ (38)

where ${\rm
gap}=\frac{\sqrt{N}}{Z}\sqrt{\sum_{n=1}^{N}\left(w_{n}^{t^{\prime}}-{\rm
E}(w^{t^{\prime}})\right)^{2}}=\frac{N}{Z}{\rm Var}(w^{t^{\prime}})=\frac{{\rm
Var}(w^{t^{\prime}})}{{\rm E}(w^{t^{\prime}})}$, with ${\rm
E}(w^{t^{\prime}})$ and ${\rm Var}(w^{t^{\prime}})$ being the empirical mean
and variance of $w^{t}$ estimated on the $N$ calibration samples. Under the
theorem assumptions, we have
$w_{n}^{t^{\prime}}=(1-\epsilon)+\epsilon\frac{\delta(S_{n}^{t^{\prime}})}{p_{C}(S_{n}^{t^{\prime}})}$.
Then, ${\rm E}(w^{t^{\prime}})=1-\epsilon+\epsilon\frac{{\rm
E}(\delta(S^{t^{\prime}}))}{p_{C}(S^{t^{\prime}})}\geq 1-\epsilon$ and ${\rm
Var}^{2}(w^{t^{\prime}})={\rm E}(w^{t^{\prime}}-{\rm
E}(w^{t^{\prime}}))^{2}\leq\epsilon^{2}{\rm
Var}^{2}(\frac{\delta(S^{t^{\prime}})}{p_{C}(S^{t^{\prime}})})\leq\frac{\epsilon^{2}}{(\min_{S}p_{C}(S^{t^{\prime}}))^{2}}{\rm
Var}^{2}(\delta(S^{t^{\prime}}))\leq\frac{\epsilon^{2}}{(\min_{S}p_{C}(S^{t^{\prime}}))^{2}}$.
The statement follows from

$\displaystyle{\rm Prob}(Z\leq a)\geq{\rm Prob}(Z\leq b)\quad{\rm if}\quad
a\geq b$ (39)

and the Quantile Lemma on $S^{t+\Delta}_{N+1}$, which implies ${\rm
Prob}(S^{t+\Delta}_{N+1}\leq Q_{C}(1-\alpha_{*}))\geq 1-\alpha_{*}$, where
$\alpha_{*}$ obeys

$\displaystyle\alpha_{*}=\alpha+{\rm gap}\leq\alpha+\frac{\sqrt{{\rm
Var}^{2}(w^{t^{\prime}})}}{{\rm
E}(w^{t^{\prime}})}\leq\frac{\epsilon}{\min_{S}p_{C}(S)}\frac{1}{1-\epsilon}.$
(40)

$\Box$

##  Appendix C Simulators

We design two different simulators: (1) a fully-synthetic simulator; and (2) a
semi-synthetic simulator.

###  C.1 General Setting

The simulators serve as fully controllable oracles to allow us test the
performance of causal inference problems. In particular, we have the following
parameters:

  * •

$N$: the total number of training users.

  * •

$M$: the total number of testing users.

  * •

$T$: the total number of steps in a time series.

  * •

$T_{0}$: the number of steps in a time series when the intervention is the
“default” one ($D=0$).

  * •

$K$: the number of different interventions (so that
$D\in\\{0,1,\cdots,K-1\\}$).

  * •

$r$: the size of latent variable dimensions.

  * •

$z_{\text{dim}}$: the size of time invariant feature dimensions.

  * •

$I$: the maximum number of non-zero intervention levels in a time series.

Whenever possible, we set the same random seeds of $1,2,3,4,5$ to aid
reproducible of our results. For the fully-synthetic simulator, a different
seed indicates that it is a different simulator (random drawn of simulator
parameters based on this seed) and also reflect the randomness comes from
neural network parameter initialization and data splitting when training
models. For the semi-synthetic simulator, the simulator parameters $\phi$ and
$\beta$ are learned from real world data, but we will need to also draw
synthetic parameters for $\psi$. In both cases, using different random seeds
can be considered to drawing new problems where each has unique parameters.

###  C.2 Fully-Synthetic Simulator

The first simulator is purely synthetic where all parameters are randomly
sampled from some pre-defined distribution. The overall data generation
process exactly follows the model we specified in Section 2.

We aims to generate a synthetic dataset for training with the following
parameters: $r=5$, $T=20$, $T_{0}=10$, $K=5$, $N=50000$, and
$z_{\text{dim}}=5$. We generate $200000$ at first place, but only use the
first $50000$ samples for training. For testing our performance, we generate
additional non-overlapping $M=5000$ samples with the same parameters, but this
time with $T=25$. We make predictions on the last $5$ time steps based on the
first $20$.

The intervention effect parameters are sampled from the following
distributions, each with a size of $K\times r$. For $d=1,2:$

$w_{1d}\sim\mathcal{U}(1,2)$ $w_{2d}\sim\mathcal{U}(1,2)$
$w_{3d}\sim\mathcal{U}(2,3).$

For $d=3,4$:

$w_{1d}\sim\mathcal{U}(-2,-1)$ $w_{2d}\sim\mathcal{U}(-2,-1)$
$w_{3d}\sim\mathcal{U}(-1,0).$

We set $w_{1d}=0$, $w_{2d}=0$ and $w_{3d}=1$ to reflect the “idle” or
“default” intervention level $d=0$. The values are designed to symmetric along
the default intervention value. We generate random time series
$\\{D_{n}^{t}\\}$ based on the identification results in Section 2.2, with the
following rules:

  1. 1.

we have enough default-intervention points ($T_{0}=10\geq r$, see Assumption
1).

  2. 2.

for each unique unit, the same type of intervention action does not repeat (by
sampling without replacement). We additionally assume at least $I>K/2$ (in
this case, $I=3$) of the interventions show up in each unit-level time series.

  3. 3.

once a $d>0$ intervention happens, there must be at least two consecutive
$d=0$ interventions (see Theorem 1).

  4. 4.

in the test dataset, at prediction time, we assume an unseen interventions in
the first $20$ time series show up at exactly $21$, and then the time series
has no intervention until the end. An example of this will be training
$D=[0,0,\cdots,1,0,\cdots,0,2,0]$ of $T=20$ and testing
$D=[0,0,\cdots,4,0,\cdots,0,3,0]+[2,0,0,0,0]$ of $T=25$.

We generate a time series $\\{\psi_{n}^{t}\\}$, based on Eq. 2, using the
corresponding $w_{1d},w_{2d},w_{3d}$ values. We generate
$\beta_{n}\sim\mathcal{N}(\mu_{n},\sigma_{n})$ for each user, each with a size
of $N\times r$, as follows:

$\mu_{n}\sim\mathcal{N}(0,1)$ $\sigma_{n}\sim\exp\\{\mathcal{N}(0,1)\\}.$

We generate unit-specific covariates $\\{Z_{n}\\}$ with $z_{\text{dim}}=5$ of
size $N\times z_{\text{dim}}$, as:

$Z_{n}\sim\mathcal{N}(0,3).$

This is then used to set up an initial $X_{n}^{0}$ (the initial starting point
of the time series) with a random-parameter MLP with input dimension of
$z_{\text{dim}}$ and output dimension of $1$. Finally, to create the synthetic
time series $\\{X_{n}^{0:T}\\}$ ($0$ comes from the initial starting point),
we randomly generate a neural network consisting of a single layer long-short
term memory (LSTM) and a MLP to generate $\phi$ at each time step. We further
compose $\phi$ with a sigmoid function so that it is bounded between $[-1,1]$.
We iteratively sample the $\\{X_{n}^{t}\\}$ points with the conditional mean
from Eq 1 (based on $X_{n}^{0:t-1}$, $Z_{n}$, and $D_{n}^{1:t}$) plus a
standard Gaussian additive error term at each step.

The outcome of this simulator gives us the following data: time series
$\\{X_{n}^{0:T}\\}$, intervention sequence $\\{D_{n}^{1:T}\\}$ and covariates
$Z_{n}$.

###  C.3 Semi-Synthetic Simulator

The second simulator is semi-synthetic, in the sense that the simulator
parameters are learned based on real world data from the WSDM Cup, organised
jointly by Spotify, WSDM, and
CrowdAI444https://research.atspotify.com/datasets/. The dataset comes from
Spotify555https://open.spotify.com/, which is an online music streaming
platform with over 190 million active users interacting with a library of over
40 million tracks. The purpose of this challenge is to understand users’
behaviours based on sequential interactions with the streamed content they are
presented with and how this contributes to skip track behaviour, which is
considered as an important implicit feedback signal.

For a detailed description of the dataset, please refer to Brost et al.
[2019]. We created our simulator based on the “Test_Set.tar.gz (14G)” file.
After unzipping the test data, we randomly chose a file named
(“log_prehistory_20180918_000000000000.csv”) as our main data source. In our
case, we only use the columns named: “session_id”, “session_position”,
“session_length”, “track_id_clean”, and “skip_3”. Here, each unique
“session_id” indicates a different sequence of interactions from an unknown
user and the length of the interaction is denoted by “session_length”, where
each step is denoted by “session_position”. The column “track_id_clean”
indicates the particular track a user listen to at each position in the
process. “skip_3” is a boolean value indicating if most of the track was
played by the user. We assign the number of $1$ if it is “True” and $0$
otherwise. We use this value to create $X_{n}^{t}$ later.

To process the data, we start by looking at the existing sequence lengths and
we select the most frequent sequence length in the data (20) (the actual
corresponding session sequence is $10$), which counts for around $57.95\%$ of
the original file. After filtering, there are around $228,460$ unique sessions
and $301,783$ unique tracks in the data. To further reduce the computational
cost of constructing a very large session-track matrix, we randomly sample a
sub-population of $6000$ unique session ids, which then bring us to $26,707$
unique tracks. We build the session-track matrix where the entry value is the
number of counts a particular track been listened over a particular session
based on this sub-population. We apply singular-value decomposition
(SVD)666https://numpy.org/doc/stable/reference/generated/numpy.linalg.svd.html
to this session-track matrix to get the session embeddings and track
embeddings. For numerical stability reason, we save only the top $50$ singular
values and further normalize it using its mean and standard deviation. For
track embeddings, we take the top $10$ singular values after the normalization
and assign them into $10$ clusters, using constrained $K$-means to to
encourage equal-sized clusters777https://github.com/joshlk/k-means-
constrained. Using these clusters, we create a new categorical variable $F$
and its continuous counterparth $\phi(F)$, which respectively represent the
corresponding cluster and its centroid vector based on the value of its
associated “track_id_clean” entry. Next, we set $X_{n}^{t}$ as the cumulative
sum of the corresponding binary “skip_3” value over each unique “session_id”
with added zero-mean Gaussian noise with standard deviation $1/2$. We keep the
first $5$ singular values for the session embedding normalization and consider
them to be the session specific covariate $Z$. We save $X$, $F$, $\phi(F)$ and
$Z$ as separate matrices as the final output, which have shapes of $6000\times
10$, $6000\times 10\times 10$, $6000\times 10\times 10$ and $6000\times 5$,
respectively.

To build a simulator, we define the following two processes:

$P(F_{n}^{t}~{}|~{}F_{n}^{1:t-1},X_{n}^{1:t-1},Z_{n})=\text{softmax}(f(F_{n}^{1:t-1},X_{n}^{1:t-1},Z_{n})),$

and

$X_{n}^{t}=\phi(F_{n}^{t})^{\mathsf{T}}\beta_{n}+\epsilon_{n}^{t},$

where $\phi(F_{n}^{t})$ is the $10$-dimensional embedding of categorical
variable $F_{n}$ given by the cluster centroid, and $\epsilon_{n}^{t}\sim
N(0,\sigma_{x}^{2})$.

Function $f(\cdot,\cdot,\cdot)$ is modeled with a deep neural network
consisting of a single layer GRU and a MLP, trained on the real $F$ with cross
entropy loss (as a likelihood based model). The second equation is trained on
the real $X$ with an ordinary least square (OLS) regression with a ridge
regularisation term ($\alpha=0.01$). Error variance is estimated by setting
$\hat{\sigma}_{x}^{2}$ to be equal to the variance of the residuals. This
defines a semi-synthetic model calibrated by our pre-processing of the real
data.

We now need to define interventions synthetically, as they are not present on
the real data. We choose $K=5$ as the size of the space of possible
interventions. We will define interventions in a way to capture the notion of
changing the exposure of tracks in particular cluster $k$ for particular user
$n$ at time $t$. Since the real data does not contain information regarding
how interventions can influence the behaviour of Spotify users, similar to the
simulator setup for the fully-synthetic one, we use $5$ random seed to
indicates different possible changes to the users behaviour in the Spotify
platform.

The intervention effect parameters are sampled from the following
distributions, each with a size of $K\times r$. For $d=1,2:$

$w_{1d}\sim\mathcal{U}(1,1.5)$ $w_{2d}\sim\mathcal{U}(1,1.5)$
$w_{3d}\sim\mathcal{U}(1.5,2.0).$

For $d=3,4$:

$w_{1d}\sim\mathcal{U}(-1.5,-1)$ $w_{2d}\sim\mathcal{U}(-1.5,-1)$
$w_{3d}\sim\mathcal{U}(0,0.5).$

We set $w_{1d}=0$, $w_{2d}=0$ and $w_{3d}=1$ to reflect the “idle” or
“default” intervention level $d=0$. The values are designed to symmetric along
the default intervention value. We generate random time series
$\\{D_{n}^{t}\\}$ based on the identification results in Section 2.2 and the
rules we adopted in the full-synthetic simulator in Appendix C.2. Notice that
the above does not correspond exactly to the model class we discuss in Section
2, as the product of intervention $\psi$ in Eq. 2 is multiplied to the softmax
operation of $f(F_{n}^{1:t-1},X_{n}^{1:t-1},Z_{n})$ and also to the
$\phi(F_{n}^{t})$.

We generate a synthetic dataset for training with the following parameters:
$r=10$, $T=35$, $T_{0}=25$, $N=3000$, and $z_{\text{dim}}=5$ (note, we have
the first $5000$ for training, but only use the first $3000$ samples are used
for actual training, as we have additional experiments to check the influence
of more training data points). For testing our performance, we generate
additional non-overlapping $M=1000$ samples with the same parameters, but this
time with $T=40$ for making predictions on the last $5$ time series based on
the first $35$.

To sample data from this simulator, we take the initial state of $X_{n}^{0}$
and $F_{n}^{0}$ from the pre-processed data. We use the $Z_{n}$ vector which
comes from the real-world user embedding. Then we iterative generated the next
$F_{n}$ via categorical draw from the softmax output based on user embedding
$Z_{n}$ all previous $F_{n}$ and $X_{n}$. Once a categorical value of $F_{n}$
is drawn, we replace it with its embedding form $\phi(F_{n})$ which is the
centroid of its cluster. We use this $\phi(F_{n})$ to generate $X_{n}$ plus an
additive Gaussian noise with zero mean and homoscedastic variance learned from
data.

The outcome of this simulator gives us the following data: time series
$\\{X_{n}^{0:T}\\}$, intervention sequence $\\{D_{n}^{1:T}\\}$ and covariates
$Z_{n}$.

##  Appendix D More Results

In this section, we present more experimental results on our method. This
includes: (1) the effect of choosing $r$; (2) further examples of
visualization and empirical examples.

###  D.1 Effect on Choice of r

The effect of different $r$ is presented in Fig. 3, with $r=5$ being the
ground truth parameter which data is generated. We have also plotted the
performance of our method under $r=3$ and $r=10$. We observed that when $r$ is
smaller than the ground truth, the model tend to under-fit the data with a
slight higher MSE error. When $r$ is bigger than the actual value, the MSE
error can be further reduced but also at the same time results in very large
variance. We also show that this can be mitigated with further regularization
such as applying $L_{1}$ and $L_{2}$ terms during likelihood learning process.

The insight we drawn from this is that when apply our approach in real-world
problems, we could potentially choose a higher $r$ without much concern about
what an exact value should be. Using a large $r$ leads to higher variance but
in general performs reasonable well compare with knowing the true $r$ (recall
that here $r=5$). We also notice that applying regularization technique even
when we know the right dimension of $r$ can be beneficial, leading to more
stable performance (with lower variance), but does not change much of the mean
estimate as shown in Fig. 3. This maybe partially caused by the gradient
optimization process when training the deep neural network. Hence, we claim
that the main price to be paid for a large choice of $r$ is that this makes
data requirements stricter (such as the number of steps prior to the first
intervention) in order to guarantee identification.

Figure 3: Effect of changing $r$ for the fully-synthetic dataset.

###  D.2 Additional Figures

Figure 4: Examples of reconstruction of training data for CSI-VAE-1 model.

Figure 5: Model-based uncertainty quantification for Reconstruction (CSI-
VAE-1).

Figure 6: Demonstration of prediction for synthetic data for CSI-VAE-1.

Figure 7: Residual distribution of our problems (based on CSI-VAE-1 model). In
general, we observe a very long tail effect across model predictions.

##  Appendix E Uncertainty Quantification

For the sake of illustration, we generate another new dataset with a random
seed of $42$ using the synthetic simulator, based on the following
specifications: $r=5$, $T=20$, $T_{0}=10$, $K=5$, $N=50000$,
$z_{\text{dim}}=5$, $I=3$. We also generated two further datasets of $M=2000$,
one for calibration and one for test.

###  E.1 Plug-in Model-based Uncertainty Quantification

We first present the model based uncertainty quantification at level of
$\alpha=0.95$, in Fig. 8. Given our conditional mean and homoscedastic error
variance $\sigma^{2}$ estimates from data, we can calculate predictive
intervals as $(\mu-1.96\times\sigma,\mu+1.96\times\sigma)$. For purely model
based uncertainty quantification with plug-in parameter estimates, we do not
need an additional calibration dataset, but it is an obviously problematic one
as it does not take into account estimation error and sampling variability.
The coverage rate is $57.75\%$ on the test dataset without further
calibration.

Figure 8: Plug-in model-based predictive interval and coverage.

###  E.2 Conformal Prediction

We then present the model-free uncertainty quantification at level of
$\alpha=0.95$ (based on the Setup 1 in section 3.2), shown in Fig. 9. We use
the absolute value over the predictive residual as our conformity score
function, $\text{abs}(y-f(x))$. We calibrate our model output based on the
calibration dataset and then apply it on the test dataset. We obtained a
coverage rate of $95.30\%$. This is a significant improvement upon what we see
in Fig. 8 and demonstrates the effectiveness of the conformal prediction
approach.

Figure 9: Conformal prediction predictive interval and coverage.