Causal Inference for Sequential Settings under Interference and Latent Confounding

Unknown

Conceptio › Archive › arXiv CS

arXiv CSopen access

Causal Inference for Sequential Settings under Interference and Latent Confounding

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

machine learning, deep learning, neural networks

Causal Inference for Sequential Settings under Interference and Latent Confounding

Phevos Paschalidis EECS, LIDS, CSAIL MIT [email protected]

Constantinos Daskalakis EECS, CSAIL MIT [email protected]

Devavrat Shah EECS, IDSS, LIDS, SDSC MIT [email protected]

arXiv:2607.14940v1 [cs.LG] 16 Jul 2026

Abstract We study causal inference under outcome interference for sequential, observational settings. Specifically, we consider settings where the binary outcomes over N units are Markovian across T time steps. At each time step, the outcomes of N units have dependencies captured through an Ising model; each outcome is also impacted through an external field capturing the effects of its treatment as well as latent confounders. Similar to panel data literature, these latent confounders are modeled to have a low-rank factor structure. Our data is a single sample from this high-dimensional distribution. To estimate causal quantities of interest, we provide a computationally efficient method based on Maximum Pseudo-Likelihood Estimation (MPLE) for learning the model parameters. Under mild assumptions, we establish non-asymptotic consistency for parameter estimation and show this translates to faithful estimation of causal quantities of interest after sampling from the learned model. We demonstrate the efficacy of the method through synthetic experiments as well as a real-world case-study investigating causal effects of vaccine rates on COVID-19 death rates within US counties nationwide.

Introduction

A key assumption underpinning much of the rich literature on causal inference is the Stable Unit Treatment Value Assumption (SUTVA), which states that the outcome of a unit is causally dependent only on its own features and treatment [1]. For many applications, however, there is interference naturally present leading to a violation of SUTVA. As a motivating example, consider a population of N individuals, some of whom receive a vaccine. Under the SUTVA assumption, each individual i’s binary outcome xi (whether or not they get sick) is dependent only on individual characteristics αi (affecting their predisposition to getting sick) and their binary treatment zi (whether they receive the vaccine). It is therefore possible to estimate the population level causal effect of the vaccine by averaging the independent outcomes of different individuals who did or did not receive it. In reality, however, each individual’s likelihood of contracting the disease is also dependent on its prevalence in the population they come in contact with. A growing body of literature has thus begun exploring causal frameworks for settings that violate SUTVA, for example [2–9]. A well-established approach is to assume that each individual’s outcome is affected by its own and others’ treatments according to an underlying network [5, 6, 9–11]. That is, xi = f (αi , zi , zN (i) ) where N (i) are the neighbors of i and f is an unknown potential outcome function (which might be randomized). An arguably more complete model of interference is that of outcome interference, which allows xi = f (αi , zi , α−i , z−i , x−i ), with the additional assumption that the joint x satisfies certain conditional independence properties capturing e.g. network dependencies [7, 8]: in the example discussed earlier, each individual’s likelihood of getting sick is directly dependent on their own predisposition, Preprint.

their own vaccination status, and the realized outcomes of their friends and family, and therefore dependent on the predispositions, vaccinations, and outcomes of the entire population. Existing work in outcome interference considers non-sequential settings [7, 8], whereas it is also useful to consider time dependence since, for example, individuals may receive the vaccine at different times and get sick at different times. In this case, we will have dependence between adjacent time steps as well. Furthermore, outcomes and interventions are impacted by confounders in observational settings, which may be latent. This is also unexplored in existing work. In summary, we have a setting where outcomes across multiple units are observed in a sequential setting, interfere with each other, and are impacted by interventions that are confounded by potentially latent variables. Contributions. We introduce a model to capture the setting of sequential outcomes over multiple units under interference and latent confounding. The latent confounding across time-steps and units is modeled to have low-rank factor structure similar to the panel data literature, cf. [12–15]. Under our model, the latent confounding manifests through the so-called unit and time specific external field parameters associated with the outcome variable. We provide a computationally efficient method based on Maximum Pseudo-Likelihood Estimation (MPLE) to learn the model parameters from a single observation of data, where recall that a single observation is a sample from a distribution of outcomes for all units and time steps with intricate dependencies. This, subsequently, enables estimation of various causal quantities of interest via sampling from the underlying distribution. For example, we can estimate the global average treatment effect which compares the average outcome over all nodes and times under complete intervention vs no intervention at all; or, more generally, we can compare the average outcomes of all nodes and times under any two choices of intervention. Synthetic data experiments and a real-world casestudy investigating the effect of vaccination on COVID-19 death rates of US counties establish the practical utility of the method. From the perspective of learning high-dimensional distributions from a single sample, our work extends a growing body of work [16–21] to the sequential setting with latent variables. From the perspective of causal inference, this work provides provable guarantees for causal estimation in a reasonably generic setting where one is given access to observational data with network interference, intertemporal dependencies, and latent confounding. While there is a rapidly growing literature on causal inference under network interference or sequential observation with latent confounding, there do not exist frameworks that offer guarantees when all complexities are co-present. 1.1

Related Works

Causal Inference under Outcome Interference. Following a surge of interest in causal inference under interference [2, 3, 22, 23], Tchetgen et al. proposed a model for outcome interference by modeling outcomes as a single draw from a Markov random field (MRF) [7]. As in our work, the MRF encodes conditional independence properties of the joint distribution of outcomes. [7] propose MPLE as a solution approach but do not provide guarantees on parameter estimation. Outcome interference √ was subsequently explored in [24, 25]. The most relevant paper to our work is [8], which explores n-consistent methods for causal effect estimation using the MRF model. Their results include guarantees on MPLE parameter estimation and mean-field inference in a non-sequential setting where confounding is from observed covariates sampled i.i.d from an underlying distribution. Causal Inference for Sequential Settings and Latent Confounding. Our approach to latent confounders is inspired by the factor model and synthetic controls literature [13–15, 26]. Unobserved confounders are unidentifiable in general [1, 27], but by considering a fixed population over time and modeling latent variables as a combination of unit-specific factors and time-specific shocks, it is possible to identify (some of the) causal quantities of interest. Under such settings, synthetic controls (and synthetic interventions) recreate the counterfactual of what would have happened had a specific unit not been treated by constructing synthetic outcomes as a linear combination of outcomes of units that did not undergo the treatment [12, 15, 28]. Recently, this line of work has been extended to account for treatment interference of the form xi = f (αi , zi , zN (i) ) [9]. All of these methods allow for counterfactual estimation of a specific type, the effect on a specific unit of always or never having the intervention, whereas our model can consider arbitrary interventional patterns. Ising Models, Exponential Family, and Latent Confounding. This work extends a line of literature on computationally efficient methods for estimating high-dimensional Ising models from a 2

single joint sample [16–21, 29–31] by allowing for sequential settings with latent variables. Prior works established parsimonious conditions under which model parameters can be learned using a single sample under generalizations of the classical Dobrushin’s uniqueness condition [32]. This collection of works builds upon use of the Maximum Pseudo-Likelihood [33, 34]. In the context of generic exponential families, computationally efficient methods for learning have been developed using different loss functions, c.f. [35–37]. This line of approach establishes its statistical properties through its connection to maximum likelihood estimation for a re-parameterized model class. This was further extended to account for unobserved confounding with a generic intervention pattern in [37]. Like this work, it connects the effect of unobserved confounding to unobserved external field. The results in [37] differ from this work in that they do not allow for interference.

Setup

We consider a network of N units exposed to binary interventions over T time steps. We denote the (t) (t) binary outcome of unit i at time t by xi ∈ {−1, 1}. Its intervention is zi ∈ {−1, 1}. In addition to its outcome and intervention, we also assume that each unit i at time t is associated with a latent (t) (t) (t) variable, αi ∈ R, which can influence both xi and zi . We are interested in a setting where our observed data is only a single realization of the trajectory (x, z) = {(x(t) , z(t) )}Tt=1 . Model. We introduce the following model to capture the temporal and spatial dependence for outcomes x(t) ∈ {−1, 1}N over t ∈ [T ] given interventions z(t) over t ∈ [T ] and latent confounders (t) αi ∈ R, i ∈ [N ], t ∈ [T ]. Precisely, we have p(x | z) =

T Y

p(x(t) |z(t) , x(t−1) ),

where

t=1

  N N N X X X X (t) (t) (t) (t) (t) (t) (t) (t−1)  , p x(t) | z(t) , x(t−1) ∝ exp  αi xi + β xi zi + ξ γij xi xj + η xi xi i=1

i=1

i̸=j

i=1

(1) for all t ∈ [T ]. We assume Γ = [γij ] ∈ RN ×N is a known symmetric matrix with zero diagonal capturing the interaction between N -dimensional outcomes at each time. Then θ = (t) (A := {αi }i∈[N ],t∈[T ] , β, ξ, η) summarize the unknown parameters, where A captures latent coun(t)

founders, β the direct effect of zi of temporal dependence.

(t)

on xi , ξ the strength of spatial dependence, and η the strength

Discussion, Key Assumptions. If ξ = 0, then (1) reduces to a more standard logistic regression model where the outcomes of each unit are independent. If ξ > 0 and Γ represents a fully connected graph, the outcome of unit i is dependent on the outcomes of all other units and, through those outcomes, the full vector of interventions and all latent features. These dependencies are propagated (t) (t) (t) along the edges of the network: xi ⊥⊥ x−i | xN (i) where N (i) = {j : γij > 0} are the neighbors of i. These conditional independencies offer statistical power even from a single sample of the joint distribution. To permit identifiability, we also need Assumptions 1-3. We write in terms of θ∗ , which denotes the true parameters of the model. (t)

Assumption 1 (Low-Rank Latent Structure). We assume that the matrix A∗ := {α∗ i }i,t ∈ RN ×T of latent confounders has rank at most k, for some k ≪ N, T . Equivalently, we can write A∗ = U ∗ (V ∗ )⊤ , where U ∗ ∈ RN ×k is a matrix of k-dimensional unit specific latent features and V ∗ ∈ RT ×k is a matrix of k-dimensional time-specific latent shocks. This low rank assumption is common in the synthetic controls and factor model literature [14, 15, 28]. Low-rank matrices have been shown to naturally arise in modern datasets and emerge from “well-behaved” generative models (e.g., Lipschitz factor functions) [38]. P Assumption 2 (Bounded Interaction). We scale ∥Γ∥∞ = maxi∈[N ] j∈[N ] |γij | = 1. We assume (t)

|ξ ∗ | < 1, maxi∈[N ],t∈[T ] |α∗ i | ≤ B, |β ∗ | ≤ B, and |η ∗ | ≤ B for some constant B ≥ 1. 3

Assumption 2 is based on Dobrushin’s uniqueness condition, which implies a number of desirable properties including concentration of measure, efficient inference, and correlation decay [39–43]. Definition 1 (Dobrushin’s Uniqueness Condition). Let µ be a measure on {−1, 1}N . Dobrushin’s interaction matrix C = (Cik )i,k∈[N ] is Cik = sup dTV µxi |x−i · | x[N ]\{i,k} , xk , µxi |x−i · | x[N ]\{i,k} , x′k . x[N ]\{i,k} xk ,x′k

Dobrushin’s coefficient is c = max1≤i≤N condition if c < 1.

k̸=i Cik . We say µ satisfies Dobrushin’s uniqueness

Specifically, the condition ∥ξ ∗ Γ∥∞ < 1 is sufficient for the conditional distribution x(t) | z(t) , x(t−1) to satisfy Definition 1 (see Lemma 2.6 of [29]). Assumption 2 also bounds the maximum influence of the latents, intervention, and temporal dependence. Given Assumptions 1 and 2, our parameter space is Θ = {(A, β, ξ, η) : A ∈ [−B, B]N ×T , |β|, |ξ|, |η| ≤ B, rank(A) ≤ k}. Assumption 3 (Excitability of Interventions). Let Z ∈ RN ×T be the matrix representation of the interventions z. Then, for any fixed θ = (A, β, ξ, η) ∈ Θ, ∥(A − A∗ ) + (β − β ∗ )Z∥2F ≥ e−cB (∥A − A∗ ∥2F + N T (β − β ∗ )2 ), for some constant c and with B from Assumption 2. In Assumption 3, we ensure that the effect of z is not completely hidden by the low-rank confounding; if Z were of sufficiently low-rank, any difference in β can be equivalently represented by a change in A which prevents identifiability. This is satisfied if, for example, z are dependent on (t) (t) confounding but still random enough such that p(zi = s|α∗ i ) ≥ e−cB for s ∈ {−1, 1}. Causal Estimands. The traditional causal estimand in the non-sequential setting under no interference is the average treatment effect, defined as ATE =

N 1 X E[xi |zi = 1] − E[xi |zi = 0] . N i=1

With the introduction of spatial and temporal dependence, the outcome of unit i at time t is dependent on the entire treatment assignment for all units at all time steps. Therefore, causal estimands of interest take the form of more general aggregated treatment effects [6]: for any two intervention patterns z0 , z1 ∈ {−1, 1}N T , define the Generalized Treatement Effect (GTE): GTE(z1 , z0 ) =

N T 1 XX (t) (t) E[xi |z = z1 ] − E[xi |z = z0 ] . N T i=1 t=1

(2)

For example, a particularly well-studied causal estimand in the context of interference is the global average treatment effect (GATE), which is effectively GTE(1, −1).

Main Results: Learning Algorithm, Guarantees

We present an algorithm for learning model parameters from single observation. The algorithm is based on Maximum Pseudo-Likelihood Estimation (MPLE). We provide non-asymptotic consistency guarantees to establish correctness by building on insights from recent works [20, 21]. Finally, we discuss its implication in terms of estimating the causal estimand defined in (2). Maximum Pseudo-Likelihood Estimation (MPLE). The maximum pseudo-likelihood is maximum likelihood estimation (MLE) applied to a conditional mean-field like approximation of true (i) (i) likelihood. Specifically, given p ≥ 1 variables Y1 , . . . , Yp and their n ≥ 1 observations y1 , . . . , yp for i ∈ [n], to estimate parameter ϕ from a parametric family of distribution over potential choices of Φ, one can use the MPLE, ϕMPLE ∈ arg max ϕ∈Φ

p n X X

(i)

log pϕ (Yj = yj |Y−j = y−j ),

i=1 j=1

where log pϕ (Yj = ·|Y−j = ·) represents the conditional distribution of Yj given all other p − 1 variables per the distribution with respect to parameter ϕ. We apply MPLE to our setting by specializing to the model of (1). The parameter of interest is θ ∈ Θ, which we estimate using the conditional distribution of x|z. Specifically, we apply MPLE (t) sequentially over t ∈ [T ], i.e. while computing the conditional likelihood of xi , we condition on (t) (t−1) (t) x−i , xi , zi , but not variables corresponding to time s for s > t. In that sense, our estimator can be viewed as sequential MPLE. Precisely, given a single sample observation x, z ∈ {−1, 1}N T , θ̂ ∈ arg min φ(θ; x, z), θ∈Θ

where

φ(θ; x, z) = −

T X N X

(t)

(t−1)

log pθ (xi |x−i , zi , xi

(3)

t=1 i=1

The benefit of MPLE in our setting is that φ is a convex function in θ and is easy to evaluate unlike MLE (since MPLE has a simple partition function). Θ is non-convex because of the low-rank restriction on A, but practical methods for low-rank optimization are well-developed [44–46]. We establish the following guarantees for the MPLE estimator. Let ∥A − A′ ∥2F + (β − β ′ )2 + (ξ − ξ ′ )2 + (η − η ′ )2 . (4) ∥θ − θ′ ∥⋆ = NT (t)

Theorem 1. Given a single observation (x, z) from (1) with parameters θ∗ = ({A∗ = [α∗ i : ˆ η̂) be the parameter estimate obtained as per (3). i ∈ [N ], t ∈ [T ]]}, β ∗ , ξ ∗ , η ∗ ), let θ̂ = (Â, β̂, ξ, Let Assumptions 1-3 hold. Then, there exists a constant C(B) = exp(O(B)) such that for any δ ∈ (0, 1/2), with probability at least 1 − δ, we have k(N + T ) log T + log 1δ . (5) ∥θ̂ − θ∗ ∥⋆ ≤ C(B) T ∥Γ∥2F We defer a full proof of Theorem 1 to Appendix B and a proof sketch to Section 5. Our dependence on 1/∥Γ∥2F is unavoidable as per Theorem 3 of [21], which extends naturally to our setting; intuitively, ∥Γ∥2F acts as our ‘effective sample size’ over the N units. Our dependence on k(N +T ) log T is from the log entropy of entry-wise bounded matrices of dimension N × T with rank k. Corollary 1 offers sufficient conditions for the error to vanish asymptotically. Corollary 1. We assume that N ≥ T and that k is a constant. If T ∥Γ∥2F = ω(N log T ), then limN,T →∞ ∥θ̂ − θ∗ ∥⋆ = 0. √ A natural setting where T ∥Γ∥2F = ω(N log T ) holds is when T = ω(1) and ∥Γ∥F = Ω( N ), e.g. if Γ is connected with bounded maximum degree. Estimating Causal Estimand. The primary causal estimand of interest is GTEθ∗ (z1 , z0 ) for any two distinct interventions z0 , z1 ∈ {−1, 1}N T . Since θ∗ is unknown, we use GTEθ̂ (z1 , z0 ) which (t) requires calculating Eθ̂ [xi |z] for all i ∈ [N ], t ∈ [T ], and z ∈ {z0 , z1 }. We use a sequential Gibbs approach which samples x(t) |z(t) , x(t−1) via Gibbs sampling for each t ∈ [T ] sequentially [47, 48]. Since the Gibbs sampler mixes fast for Ising models under Dobrushin’s condition, sequential Gibbs (t) permits efficient calculation of Eθ [xi |z] for any i ∈ [N ], t ∈ [T ], z ∈ {−1, 1}N T , and θ ∈ Θ [41]. Inference is computationally hard for models violating Dobrushin’s uniqueness condition [49]. Given estimation of GTEθ̂ (z1 , z0 ), it remains to control the pertubation error resulting from using θ̂ instead of θ∗ . Theorem 2 states a general result which leverages correlation decay for models under (t) (t) Dobrushin’s condition to argue that local errors in estimating pθ (xi = 1|z, x−i , x(−t) ) using θ′ ̸= θ do not propagate spatially and temporally [32, 43]. Assumption 2 is thus the unifying condition enabling end-to-end realization of parameter estimation, efficient inference, and pertubation bounds. Theorem 2. For any z0 , z1 ∈ {−1, 1}N T and model parameters θ and θ′ such that |η| + |ξ| < 1 and |η ′ | + |ξ ′ | < 1, there exists a constant c > 0 such that (GTEθ (z1 , z0 ) − GTEθ′ (z1 , z0 ))2 ≤ c∥θ − θ′ ∥⋆ . The proof is deferred to Appendix D. Theorems 1 and 2 imply the following. 5

ˆ < 1 and Corollary 2. Let θ̂ be the MPLE estimate of θ∗ . Under Assumptions 1-3 and if |η̂| + |ξ| ∗ ∗ |η | + |ξ | < 1, there exists a constant C(B) exponential in B such that for any δ ∈ (0, 1/2), with probability at least 1 − δ, we have k(N + T ) log T + log 1δ (GTEθ̂ (z1 , z0 ) − GTEθ∗ (z1 , z0 ))2 ≤ C(B) . T ∥Γ∥2F

Empirical Results

4.1

Synthetic Experiments

We conduct synthetic experiments to validate that our method estimates parameters and causal estimands in the presence of low-rank latent confounding with reasonable data efficiency.1 We set N = 500, T = 50, k = 3. To generate interventions, we sample W ∈ RN ×k and L ∈ RT ×k with standard Gaussian entries and define I = W ΣL⊤ , Σ = Diag(1.0, 0.7, 0.49). We normalize (t) (t) I ∈ [0, 1]N ×T and sample each zi independently with P(zi √= 1) = [I]it . We define A∗ = W Σ′ L⊤ , Σ′ = Diag(1.0, 0.8, 0.6) and normalize so that ∥A∗ ∥F / N T = 0.75. Since W and L are shared between I and A∗ , there exists latent confounding. We let Γ be Erdos-Renyi with p = 0.01, normalized so ∥Γ∥∞ = 1; we fix β ∗ = −0.3, ξ ∗ = 0.8, and η ∗ = 0.3. (0)

We sample x(0) with P(xi = 1) = 0.5 independently across i ∈ [N ]. Then, we use sequential Gibbs sampling with B = 100 local updates for each i ∈ [N ], t ∈ [T ] to generate x|z. As discussed, sequential Gibbs mixes fast since ξ ∗ < 1. We compute the MPLE parameters θ̂ by explicitly parameterizing A = U V ⊤ and minimizing L(θ, λ) = φ(θ) + λ(∥U ∥2F + ∥V ∥2F ),

(6)

where φ is from (3). This is a non-convex objective in U and V , so we alternate taking steps where U or V is fixed (see [45]). We select λ using cross-validation (see Appendix F for details) and use the same sequential Gibbs approach to estimate generalized treatment effects. Table 1 summarizes our estimation of parameters and causal estimands. We compare our approach to one that fixes ξ = 0, which we denote by θ̂ξ=0 . This is equivalent to a traditional logistic regression model that does not account for interference. We also compare to θ̂A=0 to emphasize the importance of learning latent confounding. Relative to θ̂ξ=0 and θ̂A=0 , θ̂ improves estimation of GTE(1, −1) by 92% and 91% (in absolute error), respectively. β ξ η A RMSE GTE(1, −1) θ∗ −0.300 0.800 0.300 0.000 −0.705 ± 0.010 θ̂ −0.279 ± 0.003 0.799 ± 0.019 0.293 ± 0.002 0.322 ± 0.003 −0.689 ± 0.015 θ̂ξ=0 −0.279 ± 0.003 0.000 ± 0.000 0.295 ± 0.002 0.325 ± 0.003 −0.513 ± 0.005 θ̂A=0 −0.154 ± 0.004 0.712 ± 0.026 0.222 ± 0.005 0.750 ± 0.000 −0.533 ± 0.023 Table 1: Parameter and GTE recovery. Entries report means across 10 trials with standard error. For the latent field, the reported value is the root mean square error (RMSE). We estimate GTE in each trial by averaging eight trajectories each generated with B = 100 local updates. ˆ by O((N + T ) log N T /T ∥Γ∥2 ) and Theorem 1 bounds the MSE of Â and squared error of β̂,ξ,η̂ F Theorem 2 controls the squared error of GTEθ̂ at the same rate if |ξ ∗ | + |η ∗ | < 1. Though in this case ∥Γ∥2F = 9.13 ± 0.71 and |ξ ∗ | + |η ∗ | > 1, the MPLE solution still recovers the parameters and causal effect. We also observe that latent field recovery is harder than scalar estimation, which is intuitive but not captured by our theoretical results. Nevertheless, recovering A∗ well enough to capture latent confounding and identify β permits accurate GTE estimation whereas setting A = 0 (t) fails even though α∗ i is zero in expectation. θ̂ is strictly more general than θ̂ξ=0 and, as shown in Appendix E, recovers its performance when the data satisfy no interference. 1 All

code for the synthetic and real world experiments is anonymized and available at the link https://anonymous.4open.science/r/MPLECausalInferenceC502.

4.2

COVID-19 Vaccination Case Study

We investigate the causal effects of the COVID-19 vaccine on death rates among US counties. We first validate that the observed interventional pattern permits identification. Then, we argue that our model is sufficiently expressive to capture COVID-19 outcome data. Finally, we estimate causal effects of interest and demonstrate the role of interference. Setting. We compile COVID-19 per-county death data from the NYT [50] and vaccination data from the CDC, supplemented by Bansal [51],[52]. Our outcomes are 1 if county i had more than 2 deaths per 100,000 residents during week t. Our interventions are 1 if the county was more than 30% vaccinated with a two week lag such that, e.g. vaccine rates from Sep. 1-7 affect death rates from Sep. 14-21. Thresholds are chosen to ensure diversity in x and z. The interaction graph Γ is constructed using a distance kernel function that keeps the eight nearest neighbors and weighs via ¯ exponential decay; our edge weights are wij = e−d(i,j)/d where d(i, j) is the distance between the centroids of i and j and d¯ is the median distance in the graph. Again, we normalize so ∥Γ∥∞ = 1. Our final data are the outcomes and interventions of N = 3, 014 mainland US counties with a population greater than 2, 000 from March 1st, 2020 until May 15th, 2022 (T = 115). The first county becomes vaccinated at time t = 50. Hybrid Experiments. The interventional pattern of COVID-19 vaccination may not satisfy Assumption 3 since once a county crosses the 30% vaccination threshold it remains vaccinated (see Figure E1). To validate that the data permit identifiability, we conduct hybrid experiments where z and Γ are from the data but we generate x according to our model. We do singular value decomposition on Z ∈ RN ×T and keep the features associated with the top 1 singular value (77% of energy). Then, letting k = 5, we construct W ∈ RN ×k and L ∈ RT ×k where the first columns are the features identified from Z and the other columns contain random Gaussian entries. We define A∗ = U ΣL⊤ , Σ = Diag(1.0, 0.9, 0.9, 0.7, 0.6) and normalize √ ∗ so ∥A ∥F / N T = 0.4. This represents low rank biases with partial confounding. We again fix β ∗ = −0.3, ξ ∗ = 0.8, and η ∗ = 0.3 and use sequential Gibbs with B = 100 to produce a single sample from the model. Our solution approach is unchanged from the synthetic experiments except we utilize only t ≥ 50 to calculate the gradient of β (since no county is vaccinated before t = 50). Table 2 summarizes the parameter and causal effect estimation results of θ̂, θ̂ξ=0 , and θ̂A=0 . Our results mirror the purely synthetic experiments (Table 1), demonstrating that the interventional distribution permits identification when the model is correct. β ξ η A RMSE GTE(1t≥50 , −1) θ∗ −0.300 0.800 0.300 0.000 −0.562 ± 0.001 θ̂ −0.310 ± 0.002 0.840 ± 0.001 0.299 ± 0.001 0.388 ± 0.074 −0.554 ± 0.002 θ̂ξ=0 −0.345 ± 0.017 0.000 ± 0.000 0.306 ± 0.001 0.244 ± 0.009 −0.410 ± 0.019 θ̂A=0 −0.097 ± 0.001 0.813 ± 0.004 0.278 ± 0.002 0.4 ± 0.000 −0.303 ± 0.001 Table 2: Parameter and GTE recovery for hybrid experiments. Entries report means across 10 trials with standard error. For the latent field, the reported value is the RMSE. We estimate GTE in each trial by averaging eight trajectories each generated with B = 100 local updates. Note that 1t≥50 denotes the intervention where all counties are vaccinated starting at t = 50.

Test Set Recovery. To validate that our model can capture the COVID-19 outcome data, we evaluate recovery of an unseen test set. The construction of a test set is complicated by unit and time-specific latents and dependence between units. We partition the graph and time horizon into 6 components, S6 (t) C1 , . . . , C6 and T1 , . . . , T6 , and define Ttest = j=1 {xi : i ∈ Cj , t ∈ Tj } (16.67% of data). This allows us to observe enough data from each unit and time step to reconstruct the external field (t) despite missingness. We also condition on any xi that is directly connected spatio-temporally to (r) an xj ∈ Ttest to prevent data leakage due to interference. See Figure E2 for a visualization. On the remaining outcomes (80.6% of data), we fit θ̂ and θ̂ξ=0 by minimizing (6) using cross-validation for choice of rank and λ (see Appendix F for details). After learning model parameters, we use sequential Gibbs (8 samples; B = 100) to recreate all outcomes under the observed interventions. 7

Table 3 reports the absolute error (AE) in estimating the average outcome. To offer a baseline, we repeat the same process for the hybrid setting and present analogous results in Table 4. The results indicate that the model is expressive enough to capture the real world data since the absolute error is small and on the same order of magnitude as the hybrid case. θ̂ outerperfoms θ̂ξ=0 for t ≥ 50. For θ̂ fit on the COVID data, we have ξˆ = 1.11 which violates Dobrushin’s condition. We still expect fast mixing in practice and validate by comparing results when B = 100 and B = 500 (Table E2). Train Loss Test Loss Train AE Test AE Train AE (t ≥ 50) Test AE (t ≥ 50) θ̂ 0.365 0.557 0.005 0.015 0.001 0.033 0.371 0.567 0.005 0.011 0.004 0.050 θ̂ξ=0 Table 3: Loss and absolute error for average outcome estimation for real world data. The train set comprises 279,205 outcomes (156,865 for t ≥ 50) and the test set comprises 55,226 (31,488 for t ≥ 50). The observed average outcome is −0.310 for the training set (−0.282 for t ≥ 50) and −0.312 for the testing set (−0.206 for t ≥ 50). Train Loss Test Loss Train AE Test AE Train AE (t ≥ 50) Test AE (t ≥ 50) θ̂ 0.538 0.598 0.001 0.009 0.007 0.009 θ̂ξ=0 0.551 0.613 0.001 0.011 0.006 0.019 Table 4: Loss and absolute error for average outcome estimation for hybrid data. The train and test sets are identical to the real-world setting (Table 3). The observed average outcome is 0.052 for the training set (−0.041 for t ≥ 50) and 0.053 for the testing set (−0.007 for t ≥ 50). Estimating Causal Effect of COVID-19 Vaccine. To estimate causal effects, we fit θ̂ and θ̂ξ=0 on the full data after re-selecting hyperparameters using cross validation. We use sequential Gibbs (8 samples; B = 100 local updates) to estimate counterfactuals. Figure 1 compares the observed data with the estimated evolution of the average COVID-19 outcome for two counterfactuals: a ‘no interventions’ policy where no county exceeds 30% vaccination and an ‘all interventions’ policy where all counties are vaccinated at t = 50. For t < 50, all scenarios are identical. Table 5 summarizes the estimated causal effect for t ≥ 50. For supplemental results, see Figures E3 and E4 in Appendix E.

(a) Interference.

(b) No Interference.

Figure 1: Average outcomes over time under counterfactual scenarios. As in the synthetic data, modeling interference estimates a larger effect as spill-overs enhance the efficacy of the intervention. Under interference, the ‘no interventions’ counterfactual closely follows the observed outcomes initially (when most counties were unvaccinated) before predicting worse (more positive) outcomes. The ‘all interventions’ counterfactual predicts better outcomes initially, before converging to the observed outcomes when most counties became vaccinated. Under no interference, both scenarios have similar outcomes.

Proof Sketch of Theorem 1

We divide the proof into three parts. The first two together show that, for fixed θ far from θ∗ , we have φ(θ) > φ(θ∗ ) with sufficiently high probability. The latter extends the argument to all θ ∈ Θ such that θ is far from θ∗ via a union bound and by leveraging the Lipschitzness of φ(θ). 8

‘No Interventions’ (z0 ) ‘All Interventions’ (z1 ) GTE(z1 , z0 ) θ̂ −0.189 −0.298 −0.108 θ̂ξ=0 −0.249 −0.271 −0.023 Table 5: Average outcomes for t ≥ 50 under counterfactual scenarios. Estimates average 8 samples generated using sequential Gibbs with B = 100 updates and report the standard error.

Part 1. Let θ̃ = θ − θ∗ . We rewrite our sequential MPLE as φ(θ) =

T X

φt (θ),

where

φt (θ) = −

t=1

N X

(t)

(t−1)

log pθ (xi |x−i , zi , xi

i=1

By a Taylor’s expansion, for every t ∈ [T ] there exist w ∈ [0, 1], θ′ = θ∗ + w(θ − θ∗ ) such that 1 φt (θ) − φt (θ∗ ) = θ̃⊤ ∇φt (θ∗ ) + θ̃⊤ ∇2 φt (θ′ ) θ̃. 2 ⊤ ∗ ⊤ 2 ′ Both θ̃ ∇φt (θ ) and θ̃ ∇ φt (θ ) θ̃ are random quantities dependent on x(t) |z(t) , x(t−1) which is an Ising model satisfying Dobrushin’s uniqueness condition. Therefore, for fixed t, concentration and anti-concentration results for φt (θ) − φt (θ∗ ) are known in the literature [20, 21]. We extend these results to the sequential setting with a careful martingale argument and show that T h i X E φt (θ) − φt (θ∗ ) | z(t) , x(t−1) = Ω(r(θ)), t=1 T X

h i p φt (θ) − φt (θ∗ ) − E φt (θ) − φt (θ∗ ) | z(t) , x(t−1) = O( r(θ)),

t=1

for an appropriate random function r(θ) that encodes distance of θ from θ∗ . We conclude that the undesirable event that r(θ) is large and φ(θ) − φ(θ∗ ) < 0 happens with small probability. Part 2. We show that there exists a deterministic distance notion d(θ) for which d(θ) = O(r(θ)) with high probability. In combination with Part 1, this will show that for θ with d(θ) large (i.e., θ far from θ∗ deterministically), we have φ(θ) − φ(θ∗ ) > 0 with high probability. In order to show d(θ) = O(r(θ)), we will need to argue that the “features” A, z, and x are well-conditioned with (t) (t) (t−1) high probability such that we can separate the effects of αi , βzi , and ηxi . It is at this step that we use Assumption 3 and argue that x cannot be “hidden” by A or z with high probability, which requires a concentration argument similar to Part 1. Part 3. We conclude via an ϵ-net and union bound argument that φ(θ) > φ(θ∗ ) for all θ with d(θ) sufficiently large. By the contrapositive, this implies that d(θ̂) is small which suffices to show (5).

Limitations

Our theoretical results are under Dobrushin’s uniqueness condition, which may not hold for all realworld experiments. In particular, some settings may exhibit strong outcome interference. Given our interest in causal inference, it is impossible to understand the performance of our model in estimating effects of unobserved interventions on real world data. In our experimental setting, we validate our approach by testing its ability to estimate unseen data under the implemented interventional policy.

Conclusion

We investigate causal inference for sequential settings with interference and latent confounding. We provide guarantees on parameter estimation via a computationally efficient learning approach of sequential Maximum Pseudo-Likelihood Estimation (MPLE). This enables estimation of a diverse set of causal quantities of interest. In addition to theoretical guarantees, we support our method with experiments on synthetic data and via a real-world case study investigating the effects of COVID-19 vaccines on death rates in US counties nationwide. Our work offers a new approach for causal inference in a complex setting, enabling a better understanding of causal relationships between treatments and outcomes. The establishment of causal relationships can inform important policy decisions. 9

References [1] G. W. Imbens and D. B. Rubin, Causal Inference for Statistics, Social, and Biomedical Sciences: An Introduction. Cambridge: Cambridge University Press, 2015. [2] M. G. Hudgens and M. E. Halloran, “Toward Causal Inference With Interference,” Journal of the American Statistical Association, vol. 103, pp. 832–842, June 2008. [3] P. M. Aronow and C. Samii, “Estimating Average Causal Effects Under General Interference, with Application to a Social Network Experiment,” The Annals of Applied Statistics, vol. 11, Dec. 2017. [4] F. Sävje, P. M. Aronow, and M. G. Hudgens, “Average treatment effects in the presence of unknown interference,” The Annals of Statistics, vol. 49, pp. 673–701, Apr. 2021. [5] M. D. Cattaneo, Y. He, Ruiqi, and Yu, “Robust Inference for the Direct Average Treatment Effect with Treatment Assignment Interference,” Feb. 2025. [6] V. Kandiros, C. Pipis, C. Daskalakis, and C. Harshaw, “The Conflict Graph Design: Estimating Causal Effects under Arbitrary Neighborhood Interference,” Jan. 2025. [7] E. J. Tchetgen Tchetgen, I. R. Fulcher, and I. Shpitser, “Auto-g-computation of causal effects on a network,” Journal of the American Statistical Association, vol. 116, no. 534, pp. 833–844, 2021. [8] S. Bhattacharya and S. Sen, “Causal effect estimation under network interference with meanfield methods,” The Annals of Statistics, vol. 53, no. 6, pp. 2430–2461, 2025. [9] A. Agarwal, S. H. Cen, D. Shah, and C. L. Yu, “Network Synthetic Interventions: A Causal Framework for Panel Data Under Network Interference,” Oct. 2023. [10] R. Jagadeesan, N. S. Pillai, and A. Volfovsky, “Designs for estimating the treatment effect in networks with interference,” The Annals of Statistics, vol. 48, pp. 679–712, Apr. 2020. [11] J. Ugander and H. Yin, “Randomized graph cluster randomization,” Journal of Causal Inference, vol. 11, Jan. 2023. [12] A. Abadie and J. Gardeazabal, “The economic costs of conflict: A case study of the basque country,” American economic review, vol. 93, no. 1, pp. 113–132, 2003. [13] M. Amjad, D. Shah, and D. Shen, “Robust synthetic control,” Journal of Machine Learning Research, vol. 19, no. 22, pp. 1–51, 2018. [14] A. Abadie, A. Agarwal, R. Dwivedi, and A. Shah, “Doubly robust inference in causal latent factor models,” arXiv preprint arXiv:2402.11652, 2024. [15] A. Agarwal, D. Shah, and D. Shen, “Synthetic interventions: Extending synthetic controls to multiple treatments,” Operations Research, vol. 74, no. 2, pp. 840–859, 2026. [16] S. Chatterjee, “Estimation in spin glasses: A first step,” The Annals of Statistics, pp. 1931– 1946, 2007. [17] B. B. Bhattacharya and S. Mukherjee, “Inference in ising models,” BERNOULLI, pp. 493–525, 2018. [18] C. Daskalakis, N. Dikkala, and I. Panageas, “Regression from dependent observations,” in Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, pp. 881– 889, 2019. [19] P. Ghosal and S. Mukherjee, “Joint estimation of parameters in ising model,” The Annals of Statistics, vol. 48, no. 2, pp. 785–810, 2020. [20] Y. Dagan, C. Daskalakis, N. Dikkala, and A. V. Kandiros, “Learning ising models from one or multiple samples,” in Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing, pp. 161–168, 2021. [21] V. Kandiros, Y. Dagan, N. Dikkala, S. Goel, and C. Daskalakis, “Statistical Estimation from Dependent Data,” in Proceedings of the 38th International Conference on Machine Learning, pp. 5269–5278, PMLR, July 2021. [22] C. F. Manski, “Identification of treatment response with social interactions,” The Econometrics Journal, vol. 16, no. 1, pp. S1–S23, 2013. 10

[23] M. J. van der Laan, “Causal Inference for a Population of Causally Connected Units,” Journal of causal inference, vol. 2, pp. 13–74, Mar. 2014. [24] E. Sherman and I. Shpitser, “Identification and estimation of causal effects from dependent data,” Advances in neural information processing systems, vol. 31, 2018. [25] R. Bhattacharya, D. Malinsky, and I. Shpitser, “Causal inference under interference and network uncertainty,” in Proceedings of The 35th Uncertainty in Artificial Intelligence Conference (R. P. Adams and V. Gogate, eds.), vol. 115 of Proceedings of Machine Learning Research, pp. 1028–1038, PMLR, 22–25 Jul 2020. [26] B. García Bulle, D. Shen, D. Shah, and A. E. Hosoi, “Public health implications of opening National Football League stadiums during the COVID-19 pandemic,” Proceedings of the National Academy of Sciences, vol. 119, p. e2114226119, Apr. 2022. [27] P. Ding, A first course in causal inference. Chapman and Hall/CRC, 2024. [28] A. Abadie, A. Diamond, and J. Hainmueller, “Synthetic control methods for comparative case studies: Estimating the effect of california’s tobacco control program,” Journal of the American statistical Association, vol. 105, no. 490, pp. 493–505, 2010. [29] Y. Dagan, C. Daskalakis, N. Dikkala, and S. Jayanti, “Learning from weakly dependent data under Dobrushin’s condition,” in Conference on Learning Theory, 2019. [30] C. Daskalakis, N. Dikkala, and I. Panageas, “Logistic regression with peer-group effects via inference in higher-order Ising models,” in Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, pp. 3653–3663, PMLR, June 2020. [31] C. Daskalakis, V. Kandiros, and R. Yao, “Estimating Ising Models in Total Variation Distance,” Nov. 2025. [32] P. Dobruschin, “The description of a random field by means of conditional probabilities and conditions of its regularity,” Theory of Probability & Its Applications, vol. 13, no. 2, pp. 197– 224, 1968. [33] J. Besag, “Spatial interaction and the statistical analysis of lattice systems,” Journal of the Royal Statistical Society: Series B (Methodological), vol. 36, no. 2, pp. 192–225, 1974. [34] J. Besag, “Statistical analysis of non-lattice data,” Journal of the Royal Statistical Society: Series D (The Statistician), vol. 24, no. 3, pp. 179–195, 1975. [35] A. Shah, D. Shah, and G. Wornell, “On learning continuous pairwise markov random fields,” in International conference on artificial intelligence and statistics, pp. 1153–1161, PMLR, 2021. [36] A. Shah, D. Shah, and G. Wornell, “A computationally efficient method for learning exponential family distributions,” Advances in neural information processing systems, vol. 34, pp. 15841–15854, 2021. [37] A. Shah, R. Dwivedi, D. Shah, and G. W. Wornell, “On counterfactual inference with unobserved confounding,” Sept. 2023. [38] M. Udell and A. Townsend, “Why are big data matrices approximately low rank?,” SIAM Journal on Mathematics of Data Science, vol. 1, no. 1, pp. 144–160, 2019. [39] S. Chatterjee, Concentration inequalities with exchangeable pairs. Stanford University, 2005. [40] C. Daskalakis, N. Dikkala, and G. Kamath, “Concentration of multilinear functions of the ising model with applications to network data,” Advances in Neural Information Processing Systems, vol. 30, 2017. [41] T. P. Hayes, “A simple condition implying rapid mixing of single-site dynamics on spin systems,” in Proceedings of the 47th Annual IEEE Symposium on Foundations of Computer Science, pp. 39–46, IEEE, 2006. [42] H. Künsch, “Decay of correlations under dobrushin’s uniqueness condition and its applications,” Communications in Mathematical Physics, vol. 84, no. 2, pp. 207–222, 1982. [43] H. Föllmer, “A covariance estimate for Gibbs measures,” Journal of Functional Analysis, vol. 46, pp. 387–395, May 1982. [44] B. Recht, M. Fazel, and P. A. Parrilo, “Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization,” SIAM review, vol. 52, no. 3, pp. 471–501, 2010. 11

[45] P. Jain, P. Netrapalli, and S. Sanghavi, “Low-rank matrix completion using alternating minimization,” in Proceedings of the forty-fifth annual ACM symposium on Theory of computing, pp. 665–674, 2013. [46] Y. Chi, Y. M. Lu, and Y. Chen, “Nonconvex optimization meets low-rank matrix factorization: An overview,” IEEE Transactions on Signal Processing, vol. 67, no. 20, pp. 5239–5269, 2019. [47] S. Geman and D. Geman, “Stochastic relaxation, gibbs distributions, and the bayesian restoration of images,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. PAMI6, no. 6, pp. 721–741, 1984. [48] M. J. Wainwright and M. I. Jordan, “Graphical Models, Exponential Families, and Variational Inference,” Foundations and Trends in Machine Learning, vol. 1, pp. 1–305, Dec. 2008. [49] E. Mossel, D. Weitz, and N. Wormald, “On the hardness of sampling independent sets beyond the tree threshold,” Probability Theory and Related Fields, vol. 143, no. 3, pp. 401–439, 2009. [50] The New York Times, “Coronavirus (COVID-19) data in the United States.” https:// github.com/nytimes/covid-19-data, 2021. Accessed: 2024. [51] Centers for Disease Control and Prevention, “COVID data tracker: Vaccinations.” https:// data.cdc.gov, 2021. Data endpoint: resource/8xkx-amqh. [52] S. Bansal et al., “Vaccine tracking data repository.” https://github.com/bansallab/ vaccinetracking, 2021. County-level time series vaccination data. [53] E. J. Candes and Y. Plan, “Tight oracle inequalities for low-rank matrix recovery from a minimal number of noisy random measurements,” IEEE Transactions on Information Theory, vol. 57, no. 4, pp. 2342–2359, 2011. [54] G. W. Brier, “Verification of forecasts expressed in terms of probability,” Monthly Weather Review, vol. 78, no. 1, pp. 1–3, 1950.

Representations of Causal Model

Figure A1 represents the latent confounding for a single i and t. Our consideration of latents vari(t) (t) ables in (1) is fully general for binary zi and xi . Figure A2 offers a graphical representation for the distribution over two units and time steps. (t)

αi

(t)

Figure A1: Our causal model for a single unit i at time t.

Proof of Theorem 1

In this section, we prove Theorem 1, restated below. Recall the definition of ∥ · ∥⋆ from (4). (t)

(t)

(t−1)

αi

(t−1)

(t)

(t−1)

(t)

(t−1)

(t)

αj

Figure A2: Our causal model for two units across two time steps, coupled spatially and temporally. The proof is divided into three steps: a curvature result which states that for any fixed θ sufficiently far from θ∗ , we have φ(θ) > φ(θ∗ ) with sufficiently high probability; an identifiability statement proving that we can separate the effects of the external field A from the interventions and the temporal couplings; and a union bound argument that extends the results of the first steps to all θ ∈ Θ. We outline each part in Sections B.1, B.2, and B.3. First, we define a few key terms and formalize the goals. Let θ̃ = θ − θ∗ . We also define the N -dimensional vector (t)

h(t) = [hi ]i∈[N ] ,

(t)

(t−1)

hi = αi + βzi + ηxi (t)

(t)

(7)

(t−1)

which serves as an external field for the conditional distribution x |z , x . We similarly define (t) ∗ (t) h̃ and h . We also define two notions of ‘distance’: ! ∥Ã∥2F 2 2 2 2 d(θ) = ∥θ̃∥⋆ · T ∥Γ∥F = + β̃ + ξ˜ + η̃ T ∥Γ∥2F ; (8) NT r(θ) =

T X

h i 2 ˜ (t) + h̃(t) x(t−1) , z(t) ˜ 2. E ξΓx + T ∥ξΓ∥ F 2

t=1

(9)

Note that d(θ) is deterministic and dependent only on the parameters of our model, whereas r(θ) is a random variable dependent on the data. It is relatively easy to see that a bound d(θ̂) ≤ C(B) (k(N + T ) log T + log 1/δ) immediately implies Theorem 1. Indeed, this will be our goal throughout the rest of the proof. To do so, it is helpful to first consider r(θ). In particular, the first step in our proof is to bound the undesireable event that r(θ) is large and φ(θ) fails to separate from φ(θ∗ ). Lemma 1 (Random Pointwise Bound). Let θ ∈ Θ be fixed. Then, there exist constants c1 , C1 > 0 such that for any fixed u1 > 0, we have P φ(θ) < φ(θ∗ ) + e−c1 B u21 ∩ {u21 ≤ r(θ)} ≤ C1 log N T exp(−e−c1 B u21 ); i.e. the probability that both r(θ) ≥ u21 and φ(θ) < φ(θ∗ ) + e−c1 B u21 hold is at most C1 log N T exp(−e−c1 B u21 ). 13

We dedicate Section B.1 to this result. In Section B.2, we prove r(θ) ≥ e−cB d(θ) probabilistically and then use Lemma 1 with u21 = e−cB d(θ) to show that large d(θ) implies high-probability separation of φ(θ) from φ(θ∗ ). Lemma 2 formalizes the result. Lemma 2 (Deterministic Pointwise Bound). Let θ ∈ Θ be fixed. Then, there exist constants c2 , C2 > 0 such that for any fixed u2 > 0 such that u22 ≤ d(θ), we have P φ(θ) ≥ φ(θ∗ ) + e−c2 B u22 ≥ 1 − C2 log N T exp(− e−c2 B u22 ).

In Section B.3, we show Lemma 3. Letting R denote the RHS of (10), the result follows by defining an ϵ-net over Θ = {θ = (A, β, ξ, η) : A ∈ [−B, B]N ×T , |β|, |ξ|, |η| ≤ B, rank(A) ≤ k} and proving that φ(θϵ ) > φ(θ∗ ) for all θϵ in the net such that d(θϵ ) ≥ R. Then, we extend the bound to all θ with d(θ) ≥ R using the Lipschitzness of φ(θ) and appropriate choice of ϵ in the net. Theorem 1 is a direct consequence of Lemma 3. Lemma 3 (Bound on d(θ̂)). Let θ̂ be the MPLE estimate of the true parameters. Then, there exists a constant C(B) exponentially dependent on B such that, with probability 1 − δ, we have 1 d(θ̂) ≤ C(B) k(N + T ) log T + log . (10) δ

Proof of Theorem 1. By Lemma 3, we have 1 d(θ̂) ≤ C(B) k(N + T ) log T + log , δ with probability at least 1 − δ. Recalling the definition of d(θ) in (8) and ∥ · ∥⋆ in (4), we have 1 ∥θ̂ − θ∗ ∥ · T ∥Γ∥2F ≤ C(B) k(N + T ) log T + log , δ on the same event. Dividing both sides by T ∥Γ∥2F gives the desired result. B.1

Pointwise Random Separation

We prove Lemma 1, restated below. Lemma 1 (Random Pointwise Bound). Let θ ∈ Θ be fixed. Then, there exist constants c1 , C1 > 0 such that for any fixed u1 > 0, we have P φ(θ) < φ(θ∗ ) + e−c1 B u21 ∩ {u21 ≤ r(θ)} ≤ C1 log N T exp(−e−c1 B u21 ); i.e., the probability that both r(θ) ≥ u21 and φ(θ) < φ(θ∗ ) + e−c1 B u21 hold is at most C1 log N T exp(−e−c1 B u21 ). We rewrite our sequential MPLE as φ(θ) =

T X

φt (θ),

where

φt (θ) = −

t=1

N X

(t)

(t−1)

log pθ (xi |x−i , zi , xi

i=1

By a Taylor expansion, for each fixed θ and t ∈ [T ] there exist w ∈ [0, 1], θ′ = θ∗ + w(θ − θ∗ ) such that 1 φt (θ) − φt (θ∗ ) = θ̃⊤ ∇φt (θ∗ ) + θ̃⊤ ∇2 φt (θ′ ) θ̃. (11) 2 Recalling the definition of h in (7), we write after some calculations ⊤

∗

θ̃ ∇φt (θ ) =

N X

(t)

tanh(ξ ∗ Γi x(t) + h∗ i ) − xi

i=1

(t)

˜ i x(t) ) = St , (h̃i + ξΓ

(12)

and, for ∥v∥∞ representing the maximum absolute entry for a vector v and for some constant c, N 1 ⊤ 2 1X (t) (t) ˜ i x(t) )2 θ̃ ∇ φt (θ′ ) θ̃ = sech2 (ξ ′ Γi x(t) + h′ i )(h̃i + ξΓ 2 2 i=1 N X 1 (t) 2 ′ (t) ′ (t) ˜ i x(t) )2 (h̃i + ξΓ ≥ sech (∥ξ Γx ∥∞ + ∥h ∥∞ ) 2 i=1

˜ (t) + h̃(t) ≥ e−cB ξΓx

2 2

= Ht ,

(13)

(t)

where the last inequality follows since ∥ξ ′ Γx(t) ∥∞ and ∥h′ ∥∞ are O(B). We want to argue that PT (1) , . . . , x(t) , z) and t=1 St + Ht is large with high probability. We define Gt = σ(x i 2 h ˜ 2, ˜ (t) + h̃(t) x(t−1) , z(t) + ∥ξΓ∥ (14) Q2t = E ξΓx F 2

2 t=1 Qt

which is Gt−1 -measurable; note = r(θ) where r(θ) is defined in (9). The conditional distributions of St and Ht are functions of x(t) |Gt−1 which is an Ising model satisfying Dobrushin’s uniqueness condition. Previous work has established that the conditional mean E[St + Ht |Gt−1 ] is of order Q2t and that St and Ht concentrate around their conditional means with radius Qt . To prove PT Lemma 1, we extend these arguments to the sequential setting. Lemma 4 argues that t=1 E[St + p PT P Ht |Gt−1 ] ≥ t=1 Q2t = r(θ); Lemma 5 shows that t (St + Ht ) concentrates at radius r(θ). The proofs are deferred to Appendix C. The proof of Lemma 1 synthesizes these results. Lemma 4. Let θ ∈ Θ be fixed. There exists a constant c4 > 0 such that T X

E[St + Ht |Gt−1 ] ≥ e−c4 B · r(θ).

t=1

Lemma 5. Let θ ∈ Θ be fixed. There exist constants c5 , c′5 , C5 > 0 such that for all u5 ≥ 0, ! T n p o X ′ 2 c5 B 2 P St + Ht − E[St + Ht |Gt−1 ] ≤ e max u5 r(θ), u5 ≥ 1 − C5 log N T e−c5 u5 . t=1

Proof of Lemma 1. Given u1 > 0, let Eφ = φ(θ) − φ(θ∗ ) ≥ e−c1 B u21 and Er = {u21 ≤ r(θ)}. We aim to show that P Eφc ∩ Er ≤ C1 log N T exp(−e−c1 B u21 ). As in (11), we write ∗

φ(θ) − φ(θ ) =

T X

∗

φt (θ) − φt (θ ) =

t=1

≥

T X t=1

St + Ht =

T X

1 θ̃ ∇φt (θ ) + θ̃⊤ ∇2 φt (θ′ ) θ̃ 2 ⊤

∗

E[St + Ht |Gt−1 ] +

t=1

T X

St + Ht − E[St + Ht |Gt−1 ] .

t=1 ′

By Lemmas 4 and 5, for any fixed u5 > 0 with probability at least 1 − C5 log N T e−c5 u5 , we have n p o φ(θ) − φ(θ∗ ) ≥ e−c4 B · r(θ) − ec5 B max u5 r(θ), u25 . (15) Choose u5 = e−(c4 +c5 )B u1 /2 and denote by E4∩5 the event of (15) with this choice of u5 . We claim Er ∩ E4∩5 ⊆ Eφ . It suffices to prove that when u21 ≤ r(θ), we have n p o 1 ec5 B max u5 r(θ), u25 ≤ e−c4 B r(θ) 2 15

for our choice of u5 . First, ec5 B u5

r(θ) ≤

1 −c4 B p 1 e u1 r(θ) ≤ e−c4 B r(θ). 2 2

Then, since B ≥ 1, ec5 B u25 =

1 1 −(2c4 +c5 )B 2 e u1 ≤ e−c4 B r(θ). 4 4

Therefore, (15) reduces to φ(θ) − φ(θ∗ ) ≥

1 −c4 B 1 e r(θ) ≥ e−c4 B u21 ≥ e−c1 B u21 , 2 2

c where c1 is sufficiently large. Thus, Er ∩ E4∩5 ⊆ Eφ . Then, it is also the case that Er ∩ Eφc ⊆ E4∩5 . Recalling our choice of u5 in defining E4∩5 , we have ′ −(c4 +c5 )B c e c P(E4∩5 ) ≤ C5 log N T exp − 5 u21 ≤ C1 log N T exp(−e−c1 B u21 ), 4

after selecting C1 and increasing c1 if necessary. Therefore, c P(Er ∩ Eφc ) ≤ P(E4∩5 ) ≤ C1 log N T exp(−e−c1 B u21 ).

B.2

Pointwise Deterministic Separation

We now turn towards proving Lemma 2, restated below. This translates our bound on the separation between φ(θ) and φ(θ∗ ) into one that more immediately implies parameter identification. Lemma 2 (Deterministic Pointwise Bound). Let θ ∈ Θ be fixed. Then, there exist constants c2 , C2 > 0 such that for any fixed u2 > 0 such that u22 ≤ d(θ), we have P φ(θ) ≥ φ(θ∗ ) + e−c2 B u22 ≥ 1 − C2 log N T exp(− e−c2 B u22 ). The result follows naturally from Lemma 1 and Lemma 6 below. Lemma 6. Let d(θ) and r(θ) be as defined in (8) and (9). Then, there exist constants c6 , C6 > 0 such that r(θ) ≥ e−c6 B d(θ), with probability at least 1 − C6 exp(−e−c6 B d(θ)). Proof of Lemma 2. Let u21 = e−c6 B d(θ). Lemma 1 states n o P φ(θ) ≥ φ(θ∗ ) + e−(c1 +c6 )B d(θ) ∪ {r(θ) < e−c6 B d(θ)} ≥ 1 − C1 log N T exp(−e−(c1 +c6 )B d(θ)).

(16)

By a union bound, P φ(θ) ≥ φ(θ∗ ) + e−(c1 +c6 )B d(θ) ≥ n o P φ(θ) ≥ φ(θ∗ ) + e−(c1 +c6 )B d(θ) ∪ {r(θ) < e−c6 B d(θ)} − P r(θ) < e−c6 B d(θ) . Substituting in bounds from (16) and Lemma 6 gives P φ(θ) ≥ φ(θ∗ ) + e−(c1 +c6 )B d(θ) ≥ 1 − (C1 + C6 ) log N T exp(−e−(c1 +c6 )B d(θ)). Letting c2 , C2 > 0 be suitably large and dependent on c1 , C1 , c6 , C6 gives the result since replacing d(θ) by u22 ≤ d(θ) only weakens the lower bound and the probability statement. 16

We dedicate the rest of the section to proving Lemma 6. Recalling (8), define dh and dξ as ! ∥Ã∥2F 2 2 d(θ) = + β̃ + η̃ T ∥Γ∥2F + T ξ˜2 ∥Γ∥2F . | {z } TN dξ | {z } dh

with r(θ) already defined as r(θ) =

T X

i 2 h ˜ 2. ˜ (t) + h̃(t) x(t−1) , z(t) + T ∥ξΓ∥ E ξΓx F 2

t=1

To argue that d(θ), which decouples the effects of A, β, and η, is a lower bound on r(θ), we need (t) (t) (t−1) to show that the three external-field features αi , zi , and xi are well-conditioned along the observed trajectory. The result is formalized in Lemma 7 and uses similar techniques to the proof of Lemma 1. We defer the proof to Appendix C. The proof of Lemma 6 then follows by a case analysis. (t)

Lemma 7. Let θ ∈ Θ be fixed and recall A is the N × T dimensional matrix that comprises αi . Then there exist constants c7 , C7 > 0 such that, with probability at least 1 − C7 exp − e−c7 B dh , we have T X

∥h̃(t) ∥22 ≥ e−c7 B ∥Ã∥2F + T N (β̃ 2 + η̃ 2 ) = e−c7 B dh ·

t=1

N . ∥Γ∥2F

Proof of Lemma 6. Recall that our goal is to prove r(θ) ≥ e−c6 B d(θ) with probability at least 1 − C6 exp(−e−c6 B d(θ)) for suitable constants c6 , C6 > 0. We choose a constant c sufficiently large (dependent on c7 ) and divide into cases. Case 1: dξ ≥ e−cB dh . Then, we already have ′

r(θ) ≥ dξ ≥ e−c6 B d(θ), deterministically for c′6 sufficiently large by the assumptions of our case analysis. Case 2: dh ≥ ecB dξ . Choose c large enough such that e−cB ≤ 14 e−c7 B , where c7 is the constant note that for each row i, the normalization ∥Γ∥∞ = 1 implies P from Lemma 7. Furthermore, (t) |γ | ≤ 1, so |Γ x | ≤ 1 because x(t) ∈ {−1, 1}N . Hence ∥E[Γx(t) |x(t−1) , z(t) ]∥22 ≤ N for ij i j every t. Using the elementary inequality ∥a + b∥22 ≥ 21 ∥b∥22 − ∥a∥22 , on the event of Lemma 7 it is the case that r(θ) ≥

T X t=1

T T h i 2 1 X (t) 2 X ˜ (t) + h̃(t) x(t−1) , z(t) ˜ (t) |x(t−1) , z(t) ]∥2 E ξΓx ≥ ∥h̃ ∥2 − ∥E[ξΓx 2 2 t=1 2 t=1 N 1 −c7 B ≥ e dh − dξ . ∥Γ∥2F 2

By the case assumption and our choice of c, the term in the parenthesis is at least 41 e−c7 B dh . Since ∥Γ∥2F ≤ N given our assumption ∥Γ∥∞ = 1, this implies that r(θ) ≥

′′ 1 −c7 B e dh ≥ e−c6 B d(θ), 4

on the event of Lemma 7 for appropriate choice of c′′6 . Let c6 = max{c′6 , c′′6 }. The success probability of this event is at least 1 − C6 exp(−e−c6 B d(θ)) again since dh dominates d(θ) in this case. 17

B.3

Bounding d(θ̂)

In this section, we complete the proof of Theorem 1 by translating the pointwise bound of Lemma 2 into a uniform statement over a set of all parameters in Θ = {θ = (A, β, ξ, η) : A ∈ [−B, B]N ×T , |β|, |ξ|, |η| ≤ B, rank(A) ≤ k}. The only substantive difference from a standard finite-dimensional ϵ-net argument is that the matrix component A lives in an N T -dimensional ambient space. We handle this using the low-rank constraint rank(A) ≤ k, which reduces the covering entropy to order k(N + T ). This will allow us to conclude via the contrapositive (since φ(θ̂) ≤ φ(θ∗ ) by definition) that θ̂ is close to θ∗ with high probability. In the rest of this section, we focus on the proof of Lemma 3 from which the proof of Theorem 1 follows immediately. We begin by proving a Lipschitzness result. Lemma 8. Let θ1 = (A1 , β1 , ξ1 , η1 ), θ2 = (A2 , β2 , ξ2 , η2 ) ∈ Θ. Then, |φ(θ1 ) − φ(θ2 )| ≤ 2 ∥A1 − A2 ∥1 + N T |β1 − β2 | + |ξ1 − ξ2 | + |η1 − η2 | .

Proof of Lemma 8. Define f (t) = φ(θ2 + t(θ1 − θ2 )). By the mean value theorem, d f (t) = sup |⟨θ1 − θ2 , ∇φ(θ2 + t(θ1 − θ2 ))⟩| . t∈[0,1] dt t∈[0,1]

|φ(θ1 ) − φ(θ2 )| ≤ sup

(t)

From the explicit score formula (12), each derivative with respect to an entry αi has absolute value at most 2, and the derivatives with respect to β, ξ, η are each bounded by 2N T because (t) (t−1) |zi |, |xi |, |Γi x(t) | ≤ 1. The stated bound follows immediately. We next generalize our d(θ) definition in (8) and define the deterministic distance between two arbitrary parameter values: ∥A − A′ ∥2F d(θ, θ′ ) := T ∥Γ∥2F + (β − β ′ )2 + (ξ − ξ ′ )2 + (η − η ′ )2 . (17) TN p In particular, d(θ, θ∗ ) = d(θ). Since (17) is simply a weighted Euclidean distance, d(θ, θ′ ) obeys the triangle inequality. We are now ready to prove Lemma 3. Lemma 3 (Bound on d(θ̂)). Let θ̂ be the MPLE estimate of the true parameters. Recall the definition of d(θ) in (8). Then, there exists a constant C(B) exponentially dependent on B such that, with probability 1 − δ, we have 1 d(θ̂) ≤ C(B) k(N + T ) log T + log . δ Proof of Lemma 3. Let 1 R := C(B) k(N + T ) log T + log , δ where C(B) is a sufficiently large singly-exponential function of B to be chosen below. We will show that, with probability at least 1 − δ, φ(θ) > φ(θ∗ )

for every θ ∈ Θ such that d(θ) ≥ R.

Since θ̂ minimizes φ over Θ, this immediately implies d(θ̂) ≤ R. Step 1: Constructing the net. Let NA (εA ) be an εA -net, in Frobenius norm, of the set {A ∈ [−B, B]N ×T : rank(A) ≤ k}. By a standard covering-number for rank k matrices (see Lemma 3.1 of [53]), we can bound ! √ CB N T log |NA (εA )| ≤ C k(N + T ) log . εA 18

(18)

Likewise, let Ns (εs ) be an εs -net of [−B, B]3 in Euclidean norm, so that CB . log |Ns (εs )| ≤ C log εs We define the product net N (εA , εs ) := NA (εA ) × Ns (εs ) ⊂ Θ. For every θ = (A, β, ξ, η) ∈ Θ, there exists θε = (Aε , βε , ξε , ηε ) ∈ N (εA , εs ) such that ∥A − Aε ∥F ≤ εA , ∥(β, ξ, η) − (βε , ξε , ηε )∥1 ≤ εs . p It will suffice to consider ϵA = N/T and ϵs = 1/T for which (18) and (19) give log (|N (εA , εs )|) ≤ Ck(N + T ) log (CBT ) + C log(CBT ).

(19)

(20)

Step 2: Pointwise separation on the net. By Lemma 2, for each fixed θε ∈ Θ with d(θε ) ≥ R/4, φ(θε ) ≥ φ(θ∗ ) + e−c2 B d(θε ) ≥ φ(θ∗ ) + e−c2 B R/4, except on an event of probability at most C2 log N T exp(−e−c2 B R/4). Taking a union bound over N (εA , εs ), this fails for some net point only with probability at most C2 |N (εA , εs )| log N T exp(−e−c2 B R/4).

(21)

From (20), log (C2 |N (εA , εs )| log N T ) ≤ Ck(N + T ) log (CBT ) + C log(CBT ) + log log N T ≤ ecB k(N + T ) log T, for some constant c. For C(B) chosen big enough in the expression of R, 1 exp(−e−c2 B R/4) ≤ exp(−ecB k(N + T ) log T − log ). δ Then, (21) reduces to C2 |N (εA , εs )| log N T exp(−e−c2 B R/4) ≤ δ, so we have φ(θε ) ≥ φ(θ∗ ) + e−c2 B R/4 (22) for all θε ∈ Θ such that d(θϵ ) ≥ R/4 with probability at least 1 − δ. p Step 3: Passing from the net to all of Θ. Because d(·, ·) observes the triangle inequality, for any θ ∈ Θ and nearest net point θε we have p p p d(θε ) ≥ d(θ) − d(θ, θε ). p Let C(B) ≥ 4, enlarging if necessary. Then, if d(θ) ≥ R, since ϵA = N/T , ϵs = 1/T we have 2 εA R d(θ, θε ) ≤ T ∥Γ∥2F + ε2s ≤ ∥Γ∥2F ≤ N ≤ . TN 4 Thus, d(θε ) ≥ R4 and so on the good event from Step 2 (Equation (22)) it follows that φ(θε ) ≥ φ(θ∗ ) + e−c2 B R/4. By Lemma 8, |φ(θ) − φ(θε )| ≤ 2 ∥A − Aε ∥1 + N T |β − βε | + |ξ − ξε | + |η − ηε | . √ √ Let C(B) p ≥ 50ec2 B , enlarging if necessary. Using ∥A − Aε ∥1 ≤ N T ∥A − Aε ∥F ≤ N T εA and εA = N/T , εs = 1/T , we have |φ(θ) − φ(θε )| ≤ 5N ≤ e−c2 B R/8. Therefore

φ(θ) ≥ φ(θε ) − |φ(θ) − φ(θε )| ≥ φ(θ∗ ) + e−c2 B R/8 > φ(θ∗ ).

We have thus shown that every θ ∈ Θ satisfying d(θ) ≥ R has strictly larger pseudo-likelihood than θ∗ on an event of probability 1 − δ. Since θ̂ minimizes φ over Θ, we have φ(θ̂) ≤ φ(θ∗ ). Therefore, d(θ̂) ≤ R by the contrapositive. 19

Auxiliary Lemmas

We prove the auxiliary lemmas (Lemmas 4, 5, and 7) used in the proof of Theorem 1. Lemmas 4 and 5 are used to prove Lemma 1. Lemma 7 is used in the proof of Lemma 2. Since Lemmas 1 and 7 follow from similar arguments, we begin with a proof overview. C.1

Proof Overview

PT Viewed abstractly, Lemmas 1 and 7 aim to provide a lower bound on t=1 Yt , where Yt = φt (θ) − φt (θ∗ ) in Lemma 1 and Yt = ∥h̃(t) ∥22 in Lemma 7. We choose an appropriate filtration and write T X t=1

Yt =

T X t=1

E[Yt |Ft−1 ] +

T X

T T X X Yt − E[Yt |Ft−1 ] = E[Yt |Ft−1 ] + Dt ,

t=1

where Dt = Yt − E[Yt |Ft−1 ]. Since Yt | Ft−1 is in both cases a function of an Ising model under Dobrushin’s uniqueness condition, existing results provide a lower bound for E[Yt |Ft−1 ] and a concentration result for Yt around its conditional mean (equivalently, a bound on |Dt |), for each t. PT PT A lower bound on t=1 E[Yt | Ft−1 ] follows immediately. To argue that | t=1 Dt | is sufficiently small with high probability, we prove that the conditional concentration bound for Dt implies a bound on the conditional Moment Generating Function (MGF), E[eλDt | Ft−1 ]. Then, we use our PT filtration to bound the MGF of t=1 Dt , which implies a tail bound on the sum. If the bound on the conditional MGF of Dt is in terms of a random quantity (as will be the case for Lemma 1), we will need an additional peeling argument. We state the key pieces generally. We begin with the relevant anti-concentration and concentration results for Ising models under Dobrushin’s. Lemmas 11 and Lemma 12 are Lemmas 6 and 24 in [20]; Lemma 13 follows from the same proof as Lemma 8 in [20] and is also stated in the proof of Lemma 14 in [21]. Finally, Lemma 14 follows from Theorem 4.3 of [39]. Lemma 11 (Lemma 6 of [20]). Let σ be an Ising model over {−1, 1}m satisfying Dobrushin’s uniqueness condition. For any vector a ∈ Rm , Var(a⊤ σ) ≥ c∥a∥22 . Lemma 12 (Lemma 24 of [20]). Let σ be an Ising model over {−1, 1}m with interaction matrix J and external field w satisfying Dobrushin’s uniqueness condition. Let M be a symmetric real matrix of dimension m × m with zeroes on the diagonal, let b ∈ Rm be a vector and let X f (σ) = (Mi σ + bi )(σi − tanh(Ji σ + wi )). i∈[m]

Then, for any u > 0, Pr[|f (σ)| ≥ u] ≤ exp −c min

u2 u u2 , , 2 2 ∥E[M σ + b]∥2 ∥M ∥F ∥M ∥2

Lemma 13 (Lemma 8 of [20]). Let σ be an Ising model over {−1, 1}m satisfying Dobrushin’s uniqueness condition. Let M be a symmetric real matrix of m × m with zeroes on the diagonal and let b ∈ Rm be a vector. Then, for any u > 0, Pr ∥M σ + b∥22 − E ∥M σ + b∥22 ] ≥ u c u2 ≤ exp − min ,u . ∥M ∥22 ∥M ∥2F + ∥E[M σ + b]∥22 Lemma 14 (Theorem 4.3 of [39]). Let σ be an Ising model over {−1, 1}m satisfying Dobrushin’s uniqueness condition. Let a ∈ Rm be a vector and b ∈ R be a scalar. Define, X f (σ) := ai σi + b, i∈[m]

which is an affine function of σ. Then, for any u > 0, u2 P (|f (x) − E[f (x)|]| > u) ≤ 2 exp −c · ∥a∥22

As outlined, these results will be used as subroutines to prove a lower bound for E[Yt |Ft−1 ] and conPT centration for Yt conditionally on Ft−1 . A bound on | t=1 Dt | requires first converting conditional concentration into a conditional MGF bound. Lemma 15. Let D be a real-valued random variable with E[D | F] = 0. Suppose that A and b are nonnegative F-measurable quantities satisfying A ≥ b2 , and suppose that, for all u ≥ 0, 2 u u P (|D| ≥ u | F ) ≤ C exp −c min , . A b Then there are potentially different constants c, C > 0 such that for all |λ| ≤ c/b, E eλD | F ≤ exp Cλ2 A .

(23)

Alternatively, if for all u ≥ 0 we have u2 , P (|D| ≥ u | F ) ≤ C exp −c · A then (23) holds for all λ ∈ R. PT Conditional MGF bounds for each Dt imply concentration of the sum T =1 Dt as per Lemma 16. The concentration holds for many choices of λ, which can be chosen to obtain the tightest bound. Lemma 16. Let (Dt , Ft )Tt=1 where Ft ⊆ Ft+1 are nested be such that for all t = 1, . . . , T , E[eλDt | Ft−1 ] ≤ exp(Cλ2 At ), for all λ ∈ Λ, where C > 0 is constant, At is Ft−1 -measureable, and Λ is a deterministic interval containing zero. Then, for all u, v, λ > 0 where λ ∈ Λ, we have ! T T X X 2 P Dt ≥ u, At ≤ v ≤ 2 exp(−λu + Cλ2 v 2 ). (24) t=1

t=1

Finally, if the At ’s are random and a deterministic upper bound is conservative, an additional ‘peeling’ argument is necessary. qP T Lemma 17. Suppose At are random such that 0 ≤ t=1 At ≤ J for some J > 0. Suppose also that for all u, v > 0 we have ! T T X X P Dt ≥ f (u, v), At ≤ v 2 ≤ exp(−cu2 ), (25) t=1

t=1

where f (u, v) = max{uv, bu2 } for some b ≥ 0. Then,    v u T T uX X P At  ≤ (log J) exp(−cu2 ). Dt ≥ f u, t t=1

t=1

for a potentially different constant c > 0. Note that (24) can be converted into the form of (25) by selecting λ and choosing u appropriately. The proofs of Lemmas 15, 16, and 17 are deferred to Section C.3. 21

C.2

Proof of Auxiliary Lemmas

We now complete the proofs of Lemmas 4, 5, and Lemma 7. Lemma 4 is the anti-concentration result PT lower bounding t=1 E[Yt |Ft−1 ]; Lemma 5 comprises the tail bound, MGF bound, and martingale arguments. Lemma 7 comprises all the analogous arguments, since we do not split into sub-results. Recall from (12) and (13) that St =

N X

(t)

tanh(ξ ∗ Γi x(t) + h∗ i ) − xi

(t)

˜ i x(t) ) (h̃i + ξΓ

˜ (t) + h̃(t) Ht = e−cB ξΓx

and

i=1

. 2

Recall the definitions Gt = σ(x(1) , . . . , x(t) , z), i 2 h ˜ 2 ˜ (t) + h̃(t) x(t−1) , z(t) + ∥ξΓ∥ Q2t = E ξΓx F 2

from (14), and the characterization

2 t=1 Qt = r(θ) where r(θ) is defined in (9).

Lemma 4. Let θ ∈ Θ be fixed. Then, there exists a constant c4 > 0 such that T X

E[St + Ht |Gt−1 ] ≥ e−c4 B · r(θ).

t=1

Proof of Lemma 4. First, note that h i (t) (t) (t) ˜ i x(t) )|Gt−1 , x(t) E tanh(ξ ∗ Γi x(t) + h∗ i ) − xi (h̃i + ξΓ −i (t) (t) (t) ∗ (t) ∗ (t) ˜ i x(t) ) = 0, = tanh(ξ Γi x + h i ) − E[xi |Gt−1 , x−i ] (h̃i + ξΓ (t)

where the first equality follows since Γ has zero on the diagonal so Γi x(t) is x−i measureable and (t) (t) (t) the second since E[xi |Gt−1 , x−i ] = tanh(ξ ∗ Γi x(t) + h∗ i ). Thus, by the Tower law,

E[St |Gt−1 ] = 0.

(26)

Consider now, N h i X ˜ (t) + h̃(t) ∥2 | Gt−1 = ˜ i x(t) + h̃(t) )2 |Gt−1 ] E ∥ξΓx E[(ξΓ 2 i i=1

N X

(t)

˜ i x(t) + h̃ | Gt−1 ] E[ξΓ i

˜ i x(t) |Gt−1 + Var ξΓ

i=1

˜ (t) + h̃(t) | Gt−1 ] = E[ξΓx Note that x

(t)

+ 2

N X

˜ i x(t) |Gt−1 . Var ξΓ

i=1

| Gt−1 is an Ising model under Dobrushin’s. Therefore, by Lemma 11 2 ˜ (t) + h̃(t) | Gt−1 ] + ∥ξΓ∥ ˜ 2 . E[Ht |Gt−1 ] ≥ e−c4 B E[ξΓx F 2

for some constant c4 > 0. Combining with (26) and summing over time gives the desired result by definition of r(θ) in (9). Lemma 5. Let θ ∈ Θ be fixed. There exist constants c5 , c′5 , C5 > 0 such that for all u5 ≥ 0, ! T n p o X ′ 2 P St + Ht − E[St + Ht |Gt−1 ] ≤ ec5 B max u5 r(θ), u25 ≥ 1 − C5 log N T e−c5 u5 . t=1

Proof of Lemma 5. Given Gt−1 , x(t) is an Ising model under Dobrushin’s unqiueness condition. Let ˜ J = ξ ∗ Γ, b = h̃(t) , and w = h∗ (t) . Then, t be fixed. We apply Lemma 12 to St by setting M = ξΓ, recalling the definition of Qt , for some constant c > 0 2 u u , P (|St | ≥ u | Gt−1 ) ≤ exp −c min . (27) ˜ 2 Q2t ∥ξΓ∥ Similarly applying Lemma 13 to Ht , we get for some constants c, c′ > 0, 2 ′ ˜ 2 · u | Gt−1 ≤ exp −c min u , u . P |Ht − E[Ht |Gt−1 ]| ≥ e−c B ∥ξΓ∥ ˜ 2 Q2t ∥ξΓ∥

(28)

˜ F ≥ ∥ξΓ∥ ˜ 2 , we satisfy the condition Q′ ≥ b2 . Lemma 15 therefore yields Since Qt ≥ ∥ξΓ∥ t c 2 ′2 λDt |Gt−1 ] ≤ exp(Cλ Qt ), |λ| ≤ E[e b for all t ∈ [T ]. Then, Lemma 16 gives ! T T X X 2 P Dt ≥ u, Q′ t ≤ v 2 ≤ 2 exp(−λu + Cλ2 v 2 ), t=1

t=1

for all 0 < λ ≤ cb . u c 2 2 Let λ = min{ 2Cv ≤ λu/2 so 2 , 2b }. Then, Cλ v ! 2 T T X X cu λu u 2 P Dt ≥ u, , Q′ t ≤ v 2 ≤ 2 exp − ≤ 2 exp − min . 2 4Cv 2 4b t=1 t=1

Equivalently, T X

Dt ≥ max{uv, bu },

T X

! 2 Q′ t ≤ v 2

≤ 2 exp(−cu2 ).

(29)

t=1

t=1 2

for a potentially different constant c. Since Q′t ≤ e2c0 B N for all t, applying Lemma 17 to (29) with J = C ′ N T yields  v   T T  u  uX X P Dt ≥ max ut Q′ 2t , bu2  ≤ 2 log(C ′ N T ) exp(−cu2 ).   t=1

t=1

By definition of Q′t and b, this is equivalent to  P

T X t=1

v  u T  uX ′ 2 Dt ≤ ec5 B max u5 t Qt 2 , u25  ≥ 1 − C5 log N T e−c5 u5 .    

t=1

up to a renaming of constants. The final bound follows by recalling definition of Dt and the characPT terization r(θ) = t=1 Q2t . We now turn towards Lemma 7. (t) Lemma 7. Let θ ∈ Θ be fixed and recall A is the N × T dimensional matrix that comprises αi . Then there exist constants c7 , C7 > 0 such that, with probability at least 1 − C7 exp − e−c7 B dh , we have T X

∥h̃(t) ∥22 ≥ e−c7 B ∥Ã∥2F + T N (β̃ 2 + η̃ 2 ) = e−c7 B dh ·

t=1

N . ∥Γ∥2F

Proof. We write T X

∥h̃(t) ∥22 =

t=1

T X

Lt ,

Lt :=

t=1

(t)

(t−1) 2

α̃i + β̃zi + η̃xi

i∈N

′ and let Gt′ = Gt−1 = σ(x(1) , . . . , x(t−1) , z). Furthermore, write Dt′ = Lt − E[Lt |Gt−1 ]. ′ Fix t. Since the conditional distribution of x(t−1) given Gt−1 is an Ising model with local fields bounded by O(B), there exists c > 0 such that, for both signs r ∈ {−1, 1} and every i ∈ N , (t−1) ′ Pr xi = r | Gt−1 ≥ e−cB . It follows that h i X 2 (t) (t) (t−1) 2 (t) (t) ′ E α̃i + β̃zi + η̃xi Gt−1 ≥ e−cB α̃i + β̃zi + η̃r r∈{−1,1}

(t) (t) = 2e−cB (α̃i + β̃zi )2 + η̃ 2 . Therefore, by linearity of expectation and the tower law, T X ′ E[Lt | Gt−1 ] ≥ Ce−cB ∥Ã + β̃Z∥2F + N T η̃ 2 ≥ Ce−cB ∥Ã∥2F + N T (β̃ 2 + η̃ 2 ) ,

(30)

t=1

where the second inequality follows from Assumption 3 and C, c > 0 are suitable constants. (t−1)

Since xi ∈ {−1, 1}, each single-time block contribution is an affine function of the block spin vector x(t−1) : (t−1) Lt = bt + a⊤ , t x where X (t) (t) (t) (t) at,i := 2η̃(α̃i + β̃zi ), bt := (α̃i + β̃zi )2 + η̃ 2 . i∈N

Applying Lemma 14 with f (x(t−1) ) = Lt and conditioning appropriately, we obtain for every t ∈ [T ], and u > 0, u2 ′ . Pr |Dt′ | ≥ u Gt−1 ≤ 2 exp −c′ ∥at ∥22 By Lemma 15, we therefore have ′

′ E[eλDt | Gt−1 ] ≤ exp(C ′ λ2 ∥at ∥22 ). for all λ ∈ R. Thus, Lemma 16 gives ! T T X X ′ 2 2 P Dt ≥ u, ∥at ∥2 ≤ v ≤ 2 exp(−λu + Cλ2 v 2 ), t=1

t=1

for all λ > 0. Let λ = 2Cu′ v2 . Then C ′ λ2 v 2 ≤ λu/2 so ! T T 2 X X λu ′ u ′ 2 2 P Dt ≥ u, ≤ exp −c · 2 , ∥at ∥2 ≤ v ≤ exp − 2 v t=1 t=1 for a potentially different constant c′ . Equivalently, P

T X

Dt′ ≥ uv,

t=1 (t)

t=1

! ∥at ∥22 ≤ v 2

≤ exp(−c′ u2 ).

t=1

Using (r + s)2 ≤ 2r2 + 2s2 and zi T X

T X

∈ {−1, 1}, we now note that

∥at ∥22 = 4η̃ 2

T X N X

(t)

(t) 2

α̃i + β̃zi

t=1 i=1 ≤ 8η̃ ∥Ã∥2F + 8N T η̃ 2 β̃ 2 ≤ 32B 2 ∥Ã∥2F + 16N T η˜2 + 16N T β˜2 2

′′ ≤ ec B ∥Ã∥2F + N T (β̃ 2 + η̃ 2 ) , 24

(31)

since η̃ 2 , β̃ 2 ≤ 4B. Therefore, (31) implies Pr

T X

Dt′

! q ′′ B 2 2 2 c ≥u e ∥Ã∥F + N T (β̃ + η̃ ) ≤ exp −c′ u2 .

t=1

On this event, using (30), T X

q ∥h̃(t) ∥22 ≥ Ce−cB ∥Ã∥2F + N T (β̃ 2 + η̃ 2 ) − u ec′′ B ∥Ã∥2F + N T (β̃ 2 + η̃ 2 ) .

t=1 ′

By appropriate choice of u = e−c7 B T X

ec′′ B ∥Ã∥2F + N T (β̃ 2 + η̃ 2 ) , this gives

∥h̃(t) ∥22 ≥ e−c7 B ∥Ã∥2F + N T (β̃ 2 + η̃ 2 ) = e−c7 B dh ·

t=1

N , ∥Γ∥2F

for an appropriate constant c7 > 0. The probability of the event is N ≥ 1 − C7 exp −e−c7 B dh , 1 − exp −e−c7 B dh · 2 ∥Γ∥F where the inequality follows since ∥Γ∥2F ≤ N because ∥Γ∥∞ = 1. C.3

Proof of Lemmas 15, 16, and 17

We prove Lemmas 15, 16, and 17 Lemma 15. Let D be a real-valued random variable with E[D | F] = 0. Suppose that A and b are nonnegative F-measurable quantities satisfying A ≥ b2 , and suppose that, for all u ≥ 0, 2 u u , . (32) P (|D| ≥ u | F ) ≤ C exp −c min A b Then there are potentially different constants c, C > 0 such that for all |λ| ≤ c/b, E eλD | F ≤ exp Cλ2 A .

(33)

Alternatively, if for all u ≥ 0 we have u2 P (|D| ≥ u | F ) ≤ C exp −c · , A

(34)

then (33) holds for all λ. Proof of Lemma 15. It suffices to prove the result for λ ≥ 0; the case λ < 0 follows by applying the same argument to −D, which satisfies the same tail bound. Because E[D | F ] = 0, E[eλD | F] = 1 + E[eλD − 1 − λD | F]. For y ≥ 0, by the fundamental theorem of calculus, Z y Z y Z ∞ eλy − 1 − λy = λ(eλs − 1)ds ≤ λ2 seλs ds ≤ λ2 seλs 1|y|≥s ds. 0

Using similar logic for y < 0 and via Fubini’s theorem Z ∞ E[eλD | F] ≤ 1 + λ2 seλs P(|D| ≥ s | F)ds.

(35)

Using the tail bound of (32), we rewrite (35) as E[eλD | F] ≤ 1 + Cλ2 I(λ),

Z ∞ I(λ) := 0

2 s s s exp λs − c min , ds. A b

(36)

Since s2 /A = s/b at s = A/b, I(λ) ≤ I1 (λ) + I2 (λ), where Z A/b I1 (λ) = 0

s2 s exp λs − c A

Z ∞

and

I2 (λ) =

s s exp λs − c ds. b A/b

For I1 , completing the square gives s2 λ2 A c λs − c = − A 4c A Hence ′

I1 (λ) ≤ eC λ A

Z ∞

λA s− 2c

≤ C ′ λ2 A − c ′

′ 2

′

s2 . A

se−c s /A ds ≤ c′′ AeC λ A .

For I2 , let 0 ≤ λ ≤ c/(2b). Then

s s λs − c ≤ −c′ , b b

and therefore

Z ∞ I2 (λ) ≤

′

se−c s/b ds ≤ c′′ b2 ≤ c′′ A,

0 2

since A ≥ b . Combining the two estimates gives 2

E[eλD | F] ≤ 1 + cλ2 AeCλ A ,

0 ≤ λ ≤ c/2b ′

after a renaming of constants. Writing r = λ2 A ≥ 0, the inequality 1 + creCr ≤ eC r yields E[eλD | F] ≤ exp(Cλ2 A),

0 ≤ λ ≤ c/b,

following a final renaming of constants. For the second part of the Lemma statement, note that if D instead satisfies the tail bound of (34), we can replace (36) with Z ∞ s2 λD 2 ds. E[e | F] ≤ 1 + Cλ I(λ), I(λ) := s exp λs − c A 0 The same argument gives E[eλD | F] ≤ exp(Cλ2 A) for all λ ≥ 0. Lemma 16. Let (Dt , Ft )Tt=1 be such that for all t = 1, . . . , T E[eλDt | Ft−1 ] ≤ exp(Cλ2 At )

(37)

for all λ ∈ Λ where C > 0 is constant, At is Ft−1 -measureable, and Λ is a deterministic interval containing zero. Then, for all u, v, λ > 0 where λ ∈ Λ, we have ! T T X X 2 P Dt ≥ u, At ≤ v ≤ 2 exp(−λu + Cλ2 v 2 ). t=1

Proof of Lemma 16. Let Mk =

t=1

t=1 Dt and Vk =

t=1 At . For all k ∈ [T ]:

h i E[exp(λMk − Cλ2 Vk )] = E exp λMk−1 − Cλ2 Vk−1 · exp(−Cλ2 Ak ) · E[exp(λDk )|Fk−1 ] ≤ E[exp(λMk−1 − Cλ2 Vk−1 )], 26

(38)

where the first inequality follows by the Tower law and since Ak is Fk−1 measureable and the last inequality follows by (37). Thus, E[exp(λMT −Cλ2 VT )] ≤ 1 by repeated application of (38). Now, let λ ∈ Λ with λ ≥ 0. Then, for any u and v, ! T T X X 2 P Dt ≥ u, At ≤ v = P MT ≥ u, VT ≤ v 2 t=1

t=1

≤ P λMT − Cλ2 VT ≥ λu − Cλ2 v 2

= exp(−λu + Cλ2 v 2 )E[exp(λMT − Cλ2 VT )] ≤ exp(−λu + Cλ2 v 2 ). Applying the same argument for −

t=1 Dt gives the two sided result.

Lemma 17. Suppose At are random such that 0 ≤ that for all u, v > 0 we have P

T X

Dt ≥ f (u, v),

T X

qP T

t=1 At ≤ J for some J > 0. Suppose also

! At ≤ v

≤ exp(−cu2 ),

(39)

t=1

where f (u, v) = max{uv, bu2 } for some b > 0 or f (u, v) = uv. Then,   v  u T T uX X P Dt ≥ f u, t At  ≤ (log J) exp(−cu2 ). t=1

t=1

for a potentially different constant c. qP T Proof of Lemma 17. Define qj = 2j+1 and the event Ej = {qj−1 < t=1 At ≤ qj } for j = 0, 1, . . . , log J. Then, applying (39), ! T X P Dt ≥ f (u, qj ), Ej ≤ exp(−cu2 ), t=1

qP √ T since f is monotonic in v. Since on Ej , we also have t=1 At > qj / 2, changing constants yields     v u T T uX X At  , Ej  ≤ exp(−cu2 ). P Dt ≥ f u, t t=1

t=1

A union bound over j = 0, 1, . . . , log J therefore gives   v  u T T uX X P Dt ≥ f u, t At  ≤ log(J) exp(−cu2 ). t=1

Proof of Theorem 2

D.1

Preliminaries

t=1

Theorem 2 is based on the following, which leverages Dobrushin’s uniqueness condition to argue that, due to correlation decay, error does not propagate spatially and temporally. First, for fixed z ∈ {−1, 1}N T define # " T N 1 X X (t) Mz (θ) = Eθ x z , N T t=1 i=1 i as the mean outcome under z with data generated according to θ. 27

Theorem 3. Let the interventional pattern z be fixed. If θ and θ′ are such that |η| + |ξ| < 1 and |η ′ | + |ξ ′ | < 1, there exists a constant c > 0 such that 2 Mz (θ) − Mz (θ′ ) ≤ c∥θ − θ′ ∥⋆ .

Proof of Theorem 2. Given Theorem 3, the proof is a simple application of the triangle inequality. [ 1 , z0 ))2 ≤ [Mz1 (θ) − Mz0 (θ) − (Mz1 (θ′ ) − Mz0 (θ′ ))]2 (GTE(z1 , z0 ) − GTE(z ≤ (Mz1 (θ) − Mz1 (θ′ ))2 + (Mz0 (θ) − Mz0 (θ′ ))2 ≤ 4c′ ∥θ′ − θ∥⋆ = c∥θ′ − θ∥⋆ , where c′ is the constant from Theorem 3. We therefore dedicate the section to proving Theorem 3. To do so, we invoke results from the literature on weakly dependent random variables satisfying Dobrushin’s uniqueness condition. Lemma 18 (FöllmerP Covariance Estimate). Let µ be P a Gibbs measure on {−1, 1}N with Dobrushin ∞ matrix C satisfying n≥0 C n < ∞, and let D = n=0 C n . Then for any bounded measurable functions f, g : {−1, 1}N → R, N

Covµ (f, g) ≤

1 X δi (f )Dik δk (g), 4 i,k=1

where δi (f ) := sup |f (1, x−i ) − f (−1, x−i )| x−i

is the oscillation of f at site i. Lemma 19 (Föllmer Comparison Result). Let µ and ν be probability on {−1, 1}N , asP∞ measures n sume that µ is Gibbs with Dobrushin matrix C and resolvent D = n=0 C , and define Z bk := dTV µk (· | x−k ), νk (· | x−k ) ν(dx), k ∈ [N ]. Then, for every bounded measurable f : {−1, 1}N → R, Z Z N X (bD)i δi (f ). f dµ − f dν ≤ i=1

D.2

Proof of Theorem 3

z (θw ) At a high-level, the result follows by bounding supw∈[0,1] ∂M∂ϑ for all ϑ ∈ θ (e.g., ϑ = β) and where θw = θ′ + w(θ − θ′ ). We will show that X i T s N ∂Mz (θ) 1 XX h 1 (s) (t) xi , Mϑ x(t−1) , = E Cov ∂ϑ T s=1 t=1 N i=1

(t)

where Mϑ is the multiplier of ϑ in p(x(t) |x(t−1) ) (e.g., Mβ = (s)

(t) (t) i=1 xi zi ). It therefore suffices

(t)

to control the covariances of xi , xj t ≤ s. To do so, we use Lemmas 18 and 19. Proof of Theorem 3. Let θw = θ′ + w(θ − θ′ ) for w ∈ [0, 1]. By the mean value theorem, we can write |Mz (θ) − Mz (θ′ )| =

N X T X i=1 t=1

(t)

|αi − α′ i | sup

∂Mz (θw )

w∈[0,1]

(t) ∂αi

X ϑ∈(β,ξ,η)

∂Mz (θw ) , ∂ϑ w∈[0,1]

|ϑ − ϑ′ | sup

(40) 28

Note that for any w ∈ [0, 1]. we have |ηw | + |ξw | < 1 since θw is a convex combination of θ and θ′ . We now focus on bounding the derivative of the mean outcome for all parameters. We overload notation by writing θ instead of θw . Step 1: Gradient Decomposition. We begin with the scalar parameters. Let ϑ ∈ (β, ξ, η) be fixed. We suppress conditioning on z and dependence on θ, which are assumed throughout. For arbitrary parametric models, the following holds for any h(x) not explicitly dependent on ϑ: ∂ ∂ E[h(x)] = Cov h(x), log p(x) . ∂ϑ ∂ϑ For our Markovian model, the score term decomposes as log p(x) =

T X

log p(x(t) |x(t−1) ).

t=1

Thus, also plugging in the definition Mz (θ) ! T N T ∂Mz (θ) 1X 1 X (s) X ∂ (t) (t−1) x , = Cov log p(x |x ) . ∂ϑ T s=1 N i=1 i t=1 ∂ϑ (t)

Fixing s ∈ [T ], we aim to show that the covariance term is O(1) for all ϑ. Let Mϑ be the multiplier PN (t) (t) (t) of ϑ in the conditional distribution p(x(t) |x(t−1) ). For example, Mβ = i=1 xi zi . Then, ∂ (t) (t) log p(x(t) |x(t−1) ) = Mϑ − E[Mϑ |x(t−1) ], ∂ϑ is a well-known property of exponential family distributions. For a fixed s, then, ! ! T T N N X 1 X (s) 1 X (s) X ∂ (t) (t) (t−1) (t) (t−1) Cov x , log p(x |x ) = x , Mϑ − E[Mϑ |x ] . Cov N i=1 i t=1 ∂ϑ N i=1 i t=1 (41) For each t in the summand, we condition on t − 1. Since h i (t) (t) E Mϑ − E[Mϑ |x(t−1) ]|x(t−1) = 0, Equation (41) thus reduces to " ! X X i T T N N h X X 1 1 (s) (t) (t) (t−1) (s) (t) (t−1) (t−1) E Cov . = E Cov x , Mϑ −E[Mϑ |x ]x x , Mϑ x N i=1 i N i=1 i t=1 t=1 Clearly, for t > s, this expected covariance is zero since x(t) is independent of x(s) given x(t−1) . Summarizing, X i T s N ∂Mz (θ) 1 XX h 1 (s) (t) = E Cov xi , Mϑ x(t−1) . (42) ∂ϑ T s=1 t=1 N i=1 For the remaining covariances, both across time (when t < s) and within-time (when t = s), we will use Lemmas 18 and 19. Step 2: Bounding Covariances. Step 2a: Let s = t. We note our assumptions that |ξ| < 1 and ∥Γ∥∞ = 1 imply that p(x(s) |x(s−1) ) satisfies Dobrushin’s uniqueness condition. Thus, max i

N X

Dik ≤

k=1

1 = O(1). 1 − |ξ|

(43)

We further note that (t) Mβ =

N X i=1

(t) (t) xi zi ,

1X (t) Mξ = (Γx(t) )i xi , 2 i=1 29

Mη(t) =

N X i=1

(t) (t−1)

xi xi

(t)

δk (Mϑ ) ≤ 2 (44) PN (s) (s) 1 for all ϑ ∈ (β, ξ, η). Applying Lemma 18 with f = N i=1 xi (so δi (f ) = 2/N ) and g = Mϑ (s−1) thus gives, for any fixed x , X N N N 1 1 XX 1 (s) (s) Cov xi , Mϑ x(s−1) ≤ = O(1). δi (f )Dik δk (g) ≤ N i=1 4 i=1 1 − |ξ| k=1

(t)

Step 2b: Returning now to (42), let t < s. Since Mϑ is deterministic given x(t) , we can write X N N 1 h X (s) (t) i 1 (s) (t) (t−1) (t) (t−1) = Cov . Cov x , Mϑ x E xi |x , Mϑ x N i=1 i N i=1 Both (43) and (44) still hold, so we can re-use Lemma 18 given an appropriate bound on δi (ft ), with P (s) ft = E[ i xi /N |x(t) ]. We do recursively via Lemma 19. Fix i ∈ [N ] and let y, y ′ ∈ {−1, 1}N differ only at coordinate i. Let ν = p(x(t+1) |x(t) = y ′ , z).

µ = p(x(t+1) |x(t) = y, z) Then,

# " # 1 X (s) (t) 1 X (s) (t) ′ x |x = y − E x |x = y ft (y) − ft (y ) = E N i i N i i " " # # " " # # 1 X (s) (t+1) 1 X (s) (t+1) (t) (t) ′ =E E x |x x =y −E E x |x x =y N i i N i i Z Z = ft+1 dµ − ft+1 dν. "

′

Thus, by Lemma 19, |ft (y) − ft (y ′ )| ≤

N X (bD)j δj (ft+1 ), j=1

where we have already established the base case δj (fs ) = 2/N . Since µ and ν differ only the ith coordinate, we have 0 if k ̸= i bk := tanh(|η|) ≤ |η| otherwise. Since (43) holds for any p(x(t) |x(t−1) ), we have, ′

δi (ft ) ≤ max ′ |ft (y) − ft (y )| ≤ |η| y−i ,yi ,yi

N X

Dij δj (ft+1 ) ≤

j=1

|η| 2 max δj (ft+1 ) ≤ ρs−t , (45) 1 − |ξ| j N

where ρ = |η|/(1 − |ξ|) and the last inequality follows by unravelling the recursion with the established base case. Since |η| + |ξ| < 1, we have ρ < 1. Therefore, combining (45) with (43) and (44), (t) we can use Föllmer’s covariance estimate (Lemma 18) with f = ft and g = Mϑ , which gives N N N ρs−t 1 h X (s) (t) i 1 XX (t) (t−1) (t) E xi |x , Mϑ x δi (ft )Dik δk (Mϑ ) ≤ . Cov ≤ N 4 i=1 1 − |ξ| i=1 k=1

Substituting into (42) and summing the geometric series gives X i T s N T T ∂Mz (θ) 1 XX h 1 1 X X ρs−t (s) (t) (t−1) ≤ E Cov xi , Mϑ x ≤ ∂ϑ T s=1 t=1 N i=1 T t=1 s=1 1 − |ξ| ≤

1 = cM ≤ O(1). (1 − ρ)(1 − |ξ|)

Recall, this holds for all ϑ ∈ (β, ξ, η). (t)

Step 3: Repeating for αi . Turning now to the low-rank external field, we can decompose   T N X ∂Mz (θ) 1X 1 (t−1) (s) (t)  = Cov  x , xi xi (t) T s=1 N j=1 j ∂α i

similarly to (42) and by calculating ∂ (t)

(t′ )

log p(x

(t′ −1)

∂αi

( (t) (t) xi − E[xi |x(t−1) ] if t = t′ )= 0 otherwise.

(t)

Repeating Step 2 with g = xi and δj (g) = 2(I{j = i}) gives ∂Mz (θ) cM ≤ . (t) N T ∂α i

z (θ) Step 4: Completing the Argument. We have thus shown that ∂M∂ϑ = cM for ϑ ∈ (β, ξ, η) and ∂Mz (θ) = cM /N T . Returning to our expression of the difference in average treatment estimates in ∂ϑ (40) and using Cauchy-Schwartz,

|Mz (θ) − Mz (θ′ )| ≤

N T X cM X X (t) (t) |αi − α′ i | + cM |ϑ − ϑ′ | N T i=1 t=1 ϑ∈(β,ξ,η) v u T N X u 1 X (t) (t) ≤ cM t (α − α′ i )2 + cM N T i=1 t=1 i

Thus, (Mz (θ) − Mz (θ′ ))2 ≤ c

∥A − A′ ∥2F +c NT

Supplemental Experimental Results

E.1

Performance under No Interference

(ϑ − ϑ′ )2 .

ϑ∈(β,ξ,η)

(ϑ − ϑ′ )2 ≤ c∥θ − θ′ ∥⋆ .

ϑ∈(β,ξ,η)

We conduct a further synthetic experiment to investigate the performance of θ̂ when the underlying distribution satisfies no interference. Our experimental setting is identical to Section 4.1 except for fixing ξ ∗ = 0. Table E1 summarizes the parameter and causal effect estimation results. Our results show that θ̂ is strictly more general than θ̂ξ=0 since it correctly identifies the lack of interference. θ∗ θ̂ θ̂ξ=0 θ̂A=0

β −0.300 −0.280 ± 0.003 −0.280 ± 0.003 −0.155 ± 0.004

ξ 0.000 −0.001 ± 0.020 0.000 ± 0.000 −0.005 ± 0.030

η 0.300 0.296 ± 0.001 0.296 ± 0.001 0.225 ± 0.005

A RMSE 0.000 0.320 ± 0.003 0.320 ± 0.003 0.750 ± 0.000

GTE(1, −1) −0.534 ± 0.002 −0.516 ± 0.006 −0.515 ± 0.005 −0.371 ± 0.010

Table E1: Parameter and GTE recovery. Entries report means across 10 trials with standard error. For the latent field, the reported value is the root mean square error (RMSE). We estimate GTE in each trial by averaging eight trajectories each generated with B = 100 local updates.

E.2

Supplement to Real-World Case Study

Interventional Distribution. Figure E1 plots the percentage of US counties above the 30% threshold over time. The first county crosses the threshold at time t = 50 and by t = 90 more than 90% of counties observe the intervention. 31

Figure E1: The percentage of US counties above the 30% vaccination threshold over time.

E.3

Test Set Construction

Figure E2 offers a visualization of our test set construction on the grid graph as described in Section 4.2. We partition the units (y-axis) and time (x-axis) into 6 sets. Nodes in the test set are shaded black. Nodes that are conditioned on to prevent data leakage due to interference are shaded gray; these are all nodes that are directly connected spatio-temporally to a node in the test set. Nodes in the training set are white. Appendix F.1 discusses construction of cross-validation sets, which is by a similar procedure. Note that in the presence of the unseen test set, cross-validation for hyperparameter selection is performed on the training set which is further partitioned (see Appendix F.1).

Figure E2: A visualization of our procedure for constructing a held-out test set. We assume the graph is a grid graph so nodes are connected contiguously. Nodes shaded black, gray, and white comprise the test set, the conditioned upon “separator set”, and the training set, respectively. Validating MCMC Mixing. To validate that sequential Gibbs with B = 100 is sufficiently mixed for θ̂ fit on COVID-19 data, we re-compute the statistics of Table 3 using B = 500. Table E2 compares the MCMC estimate of the average outcome in the test and train sets when B = 100 and B = 500 for θ̂ and θ̂ξ=0 . Since the latter satisfies no interference, the sequential Gibbs approach is equivalent to sampling independent Bernoullis and mixes immediately. We expect and observe negligible differences between inference with B = 100 and B = 500. We observe similarly neglible differences in estimates for θ̂ at B = 100 and B = 500, suggesting that B = 100 suffices to sample from the Ising model conditional distribution. Supplement to Causal Inference. Figure E3 visualizes the estimated average COVID-19 outcome per-unit under the two counterfactuals for θ̂ and θ̂ξ=0 , similar to Figure 1. The results reinforce that modeling interference estimates a larger causal effect of the COVID-19 vaccine. 32

θ̂ θ̂ξ=0

B 100 500 100 500

Train AE 0.005 0.004 0.005 0.005

Test AE 0.015 0.013 0.011 0.011

Train AE (t ≥ 50) 0.001 0.001 0.004 0.003

Test AE (t ≥ 50) 0.033 0.034 0.050 0.047

Table E2: Absolute error for average outcome estimation for real world data when B = 100 and B = 500.

(a) Interference.

(b) No Interference.

Figure E3: Average outcomes per-unit under counterfactual scenarios as estimated by θ̂ (interference) and θ̂ξ=0 (no interference).

As a final validation, we also simulate outcomes under the observed policy and compare the estimated per-time and per-unit averages with the observed data for θ̂ and θ̂ξ=0 . Figure E4 reinforces that the model is sufficiently expressive to capture COVID-19 outcomes since the per-time and per-unit outcomes closely match the data.

Figure E4: Estimated average outcomes per-time and per-unit under the observed intervention compared to the observed data. Both θ̂ and θ̂ξ=0 closely follow the observed data.

Experimental Details

F.1

Cross-Validation

We use a cross-validation approach for selecting hyperparameters by creating validation sets similar to our construction of the test set. Figure E2 offers a visualization of our construction if we were to partition into five sets (whereas we use 6 for test set construction and 7 for cross validation). F.1.1

Set Construction (t)

Denote by T = {xi : i ∈ [N ], t ∈ [T ]} the set of all outcomes. There are two construction approaches, based on whether Ttrain = T . For the synthetic and hybrid experiments and when creating counterfactual estimates on real world data, we train on the full data. For the test set recovery experiments, Ttrain ⊂ T . 33

If Ttrain = T , we partition the nodes and time horizon into 7 subsets C1 , . . . , C7 and T1 , . . . , T7 . We choose the 7 validation sets as 7 [ (t) Vk = {xi : i ∈ Cj , t ∈ Tj+k }, k = 0, . . . , 6, j=1

with indices wrapping around. We define by Sk the nodes that are directly connected spatiotemporally with a node in the validation set. That is, n o (t) (t−1) (t+1) (t) Sk = xi ∈ / Vk : xi ∈ Vk or xi ∈ Vk or ∃j with γij > 0 and xj ∈ Vk . For fold k, we holdout both Vk and Sk and fit on Fk = Ttrain \ {Vk ∪ Sk }. Sk acts as a Markov blanket that enforces conditional independence between Vk and Fk during training and evaluation. If Ttrain ⊂ T , the construction procedure must also respect independence between the test and training sets. Let Ttest denote the test set, o n (t) (t−1) (t+1) (t) Tsep = xi ∈ / Ttest : xi ∈ Ttest or xi ∈ Ttest or ∃j with γij > 0 and xj ∈ Ttest , denote the separator set for the test set, and Ttrain = T \ {Ttest ∪ Tsep } be the training set. By construction, the nodes whose outcomes comprise the train set are not the same over time; i.e., for (t) (r) t ̸= r we may have xi ∈ Ttrain but xi ∈ / Ttrain . Therefore, we cannot create a stable partition over the nodes in Ttrain to repeat the validation construction. Instead, we create dynamic partitions of the nodes and combine over time conservatively. We denote (t) by Pt = {i ∈ [N ] : xi ∈ Ttrain } the available nodes at time t and partition each Pt into 7 sets Ct,1 , . . . , Ct,7 . If Pt = Pr , we use the same partition. In general, the test set construction is such that Pt = Pt−1 except for a few time steps when the set of nodes comprising the test set changes. We partition the time horizon into T1 , . . . , T7 and exclude the few times t for which Pt ̸= Pt−1 . Thus, Ct,1 , . . . , Ct,7 are dynamic partitions of the nodes and T1 , . . . , T7 build in a buffer when the partitions change. We define the validation sets as Vk =

7 [

(t)

{xi : i ∈ Ct,j , t ∈ Tj+k },

k = 0, . . . , 6.

j=1

Clearly, Vk ⊆ Ttrain . The separator set is o n (t) (t−1) (t+1) (t) Sk = Ttrain ∩ xi ∈ / Vk : xi ∈ Vk or xi ∈ Vk or ∃j with γij > 0 and xj ∈ Vk . Finally, we define Fk = Ttrain \ {Vk ∪ Sk }. Sk ensures the independence of Vk and Fk during training and evaluation. Tsep separates the test set from Vk ∪ Sk ∪ Fk and ensures independence from the cross-validation process. F.1.2

Evaluation

In both cases, we run an MPLE fit for each combination of hyperparameters on each validation fold. Since Sk separates Fk and Vk , the pseudo-likelihood of Fk is independent from Vk . To evaluate performance, we do conditional inference by conditioning on the observed outcomes of Sk and using sequential Gibbs (16 samples; B = 10) to estimate the average outcome of the held out validation set. Conditional on Sk , the outcomes of Vk are independent from Fk . We select the candidate with the lowest absolute error in estimating the average outcome of Vk in the post-interventional period (i.e., considering only t ≥ s where s is the first time at which a node observes the intervention). Note s = 1 for synthetic experiments and s = 50 for hybrid and real-world experiments. If multiple fits are within a standard error of the best fit, we select among them the candidate with the lowest conditional Brier score [54], defined as the squared difference between the predicted probability of each outcome and its occurrence. Formally, X 1 (t) (t) (t) 2 Brier(θ) = Pθ (xi = 1|x−i , z(t) , x(t−1) ) − xi . |Vk | (t) (i,t):xi ∈Vk

If multiple eligible candidates are within a standard error of the best conditional brier score, we choose the least regularized option. 34

F.1.3

Results

Let Λ = {0.001, 0.005, 0.01, 0.05, 0.1, 0.5}, which we use for all searches. Below, we list the hyperparameters used for all results. See the repository for full cross-validation results. • For the synthetic experiment of Section 4.1, we fix k = 3. After evaluation, we select λ = 0.05 for both θ̂ and θ̂ξ=0 . • For the hybrid experiment, we fix k = 5. After evaluation, we select λ = 0.001 for θ̂ and λ = 0.05 for θ̂ξ=0 . • For the test set recovery experiment for real-world data, we perform a grid search over k ∈ {3, 5, 8} and λ ∈ Λ. We select k = 8 and λ = 0.01 for both θ̂ and θ̂ξ=0 . • For the test set recovery experiment for hybrid data, we fix k = 5 and select λ = 0.005 for both θ̂ and θ̂ξ=0 . • For causal effect estimation on real-world data, we re-do the grid search over k ∈ {3, 5, 8} and λ ∈ Λ using T = Ttrain . We select k = 5 and λ = 0.05 for θ̂ and k = 5 and λ = 0.01 for θ̂ξ=0 . • For the synthetic results under no interference (Section E.1), we fix k = 3 and select λ = 0.05 for both θ̂ and θ̂ξ=0 . F.2

Computational Resources

All data generation, optimization, and sampling are run on a cluster with multiple CPU cores. Sampling may take several hours. Results can be reproduced on a standard desktop/laptop.

Record · ID 373395 · SHA-256 07961240465ab178
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.