# Inferring the dependence graph density of binary graphical models in high dimension

Julien Chevallier, Eva Löcherbach, Guilherme Ost

<figure><img src="image_1.png" /></figure>

## To cite this version:

Julien Chevallier, Eva Löcherbach, Guilherme Ost. Inferring the dependence graph density of binary graphical models in high dimension. 2024. hal-04609972

HAL Id: hal-04609972

https://hal.science/hal-04609972v1

Preprint submitted on 18 Apr 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
INFERRING THE DEPENDENCE GRAPH DENSITY OF BINARY
GRAPHICAL MODELS IN HIGH DIMENSION

BY JULIEN CHEVALLIER¹,ª ID, EVA LÖCHERBACH²,ᵇ ID AND GUILHERME OST³,ᶜ ID

¹Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP*, LJK, 38000 Grenoble, France,

ªjulien.chevallier1@univ-grenoble-alpes.fr

²Université Paris 1 Panthéon Sorbonne, ᵇEva.Locherbach@univ-paris1.fr

³Institute of Mathematics, Federal University of Rio de Janeiro, ᶜguilhermeost@im.ufrj.br

We consider a system of binary interacting chains describing the dynamics of a group of *N* components that, at each time unit, either send some signal to the others or remain silent otherwise. The interactions among the chains are encoded by a directed Erdös-Rényi random graph with unknown parameter *p* ∈ (0, 1). Moreover, the system is structured within two populations (excitatory chains versus inhibitory ones) which are coupled via a mean field interaction on the underlying Erdös-Rényi graph. In this paper, we address the question of inferring the connectivity parameter *p* based only on the observation of the interacting chains over *T* time units. In our main result, we show that the connectivity parameter *p* can be estimated with rate $N^{-1/2} + N^{1/2}/T + (\log(T)/T)^{1/2}$ through an easy-to-compute estimator. Our analysis relies on a precise study of the spatio-temporal decay of correlations of the interacting chains. This is done through the study of coalescing random walks defining a backward regeneration representation of the system. Interestingly, we also show that this backward regeneration representation allows us to perfectly sample the system of interacting chains (conditionally on each realization of the underlying Erdös-Rényi graph) from its stationary distribution. These probabilistic results have an interest in its own.

CONTENTS

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Model definition, notation and main results</td><td>5</td></tr><tr><td>2.1</td><td>Heuristics for the spatio-temporal mean</td><td>8</td></tr><tr><td>2.2</td><td>Heuristics for the spatial variance</td><td>10</td></tr><tr><td>2.3</td><td>Heuristics for the temporal variance</td><td>11</td></tr><tr><td>3</td><td>Backward regeneration scheme</td><td>13</td></tr><tr><td>3.1</td><td>Backward regeneration representation</td><td>13</td></tr><tr><td>3.2</td><td>Coalescence of the backward random walks</td><td>15</td></tr><tr><td>4</td><td>Key steps and proof of Theorem 2.1</td><td>16</td></tr><tr><td>5</td><td>Simulation study</td><td>18</td></tr><tr><td>5.1</td><td>Practical implementation</td><td>18</td></tr><tr><td>5.2</td><td>Choice of the tuning parameter</td><td>18</td></tr><tr><td>5.3</td><td>Numerical results</td><td>19</td></tr><tr><td></td><td>Acknowledgments</td><td>21</td></tr><tr><td></td><td>Funding</td><td>21</td></tr><tr><td></td><td>References</td><td>21</td></tr></table>

<footnote>*Institute of Engineering Univ. Grenoble Alpes</footnote>

<footnote>MSC2020 subject classifications: Primary 62M05; secondary 60J10, 60K35, 62F12.</footnote>

<footnote>Keywords and phrases: Dependence graph inference, Markov chain, mean field limit, perfect simulation.</footnote>
<table><tr><td>A</td><td>Proof of Theorem 3.3</td><td>23</td></tr><tr><td>B</td><td>Coalescence couplings</td><td>24</td></tr><tr><td>B.1</td><td>Two main lemmas</td><td>24</td></tr><tr><td>B.2</td><td>Construction of the couplings</td><td>24</td></tr><tr><td>B.3</td><td>Properties of the couplings</td><td>26</td></tr><tr><td>B.4</td><td>Proof of the two lemmas</td><td>29</td></tr><tr><td>C</td><td>Coalescence of two or more backward random walks</td><td>29</td></tr><tr><td>D</td><td>Proof of Lemma 3.8</td><td>35</td></tr><tr><td>E</td><td>Temporal convergences</td><td>38</td></tr><tr><td>E.1</td><td>Proof of Equation (32)</td><td>38</td></tr><tr><td>E.2</td><td>Proof of Equation (34)</td><td>39</td></tr><tr><td>E.3</td><td>Proof of Equation (33)</td><td>51</td></tr><tr><td>F</td><td>Proof of Proposition 4.2</td><td>54</td></tr><tr><td>F.1</td><td>General notation</td><td>55</td></tr><tr><td>F.2</td><td>Study of the rescaled random environment</td><td>55</td></tr><tr><td>F.2.1</td><td>Notation</td><td>56</td></tr><tr><td>F.2.2</td><td>Convergence rates</td><td>56</td></tr><tr><td>F.3</td><td>Results regarding the inverse matrix</td><td>57</td></tr><tr><td>F.3.1</td><td>Notation</td><td>57</td></tr><tr><td>F.3.2</td><td>Convergence rates</td><td>58</td></tr><tr><td>F.4</td><td>Proof of Lemma F.1</td><td>59</td></tr><tr><td>F.5</td><td>Proof of Lemma F.4</td><td>62</td></tr><tr><td>F.6</td><td>Proof of Lemma F.6</td><td>68</td></tr><tr><td>F.7</td><td>Proof of Proposition 4.2</td><td>70</td></tr><tr><td>G</td><td>Inversion of Ψ</td><td>74</td></tr><tr><td>H</td><td>Auxiliary results</td><td>77</td></tr><tr><td>I</td><td>A note on lower bounds</td><td>79</td></tr><tr><td>I.1</td><td>Statistical setting 1 and problem formulation</td><td>79</td></tr><tr><td>I.2</td><td>Statistical setting 2 and problem formulation</td><td>81</td></tr><tr><td>I.3</td><td>Auxiliary results</td><td>83</td></tr></table>

**1. Introduction.** Understanding how to infer and interpret interactions between measured components of large scale systems is a central question in many scientific fields, including Neuroscience, Statistical Physics, and Social Networks (Strogatz, 2001). This has been done traditionally in the framework of probabilistic graphical models (Lauritzen, 1996), in which the interaction among the components of the system is encoded through a graph, sometimes referred to as the graph of conditional dependencies. For such models, a fundamental question is that of estimating the underlying graph of conditional dependencies or some function of it from data. Over the past three decades, this question has been extensively investigated not only under the assumption that the data is a set of independent and identically distributed observations and the underlying graphical model has pairwise interactions (Bresler, 2015; Montanari and Pereira, 2009; Ravikumar, Wainwright and Lafferty, 2010), but also for time-dependent data (Eichler, 2012; Duarte et al., 2019; Reynaud-Bouret, Rivoirard and Tuleau-Malot, 2013) or graphical models with high-order interactions (Basu and Michailidis, 2015; Lerasle and Takahashi, 2016). In recent years, understating the theoretical guarantees of the proposed methods (e.g., consistency, computational complexity) in the high-dimension setting became a major issue, as now the simultaneous activity of many components can be routinely recorded. In high-dimension, most of the works assume that the underlying graph of dependencies is sparse in a suitable sense, showing that under this assumption the proposed method performs reasonably well. The cases in which the graphs of
dependencies are dense, on the other hand, are much less studied. In this setting we are only aware of the works of (Byol Kim and Kolar, 2021; Delattre and Fournier, 2016) and (Liu, 2020). In these works the focus is not on the estimation of the graph of dependencies of the corresponding graphical model, as we detail in what follows.

In (Byol Kim and Kolar, 2021), particularly motivated by applications in Genetics and Neuroscience, the authors propose a density ratio approach to estimate the “difference network”, the difference between two graphs of dependencies. The key assumption in their approach is that this difference network is sparse and, in particular, their method can also be applied when each individual graph is dense as long as their difference is sparse. (Delattre and Fournier, 2016) consider a graphical model with pairwise interactions for which the graph of dependencies is a realization of a weighted Erdős-Rényi random graph with edge probability p. Specifically, (Delattre and Fournier, 2016) work with N Hawkes point processes on $[0, T]$ - a class of point processes in which the intensity at any given time is a linear function of the past events - where each pair of these point processes are independently coupled with probability p. Moreover, the existing interactions are excitatory and of mean-field type (i.e., the occurrence of an event increases the chance of a new event to occur and the increment scales as $1/N$). The main goal of the paper is not to estimate the random graph of dependencies but rather the density of connections of this graph (the parameter p) based only on the information present in the N point processes on $[0, T]$. Under some few extra assumptions, the authors show that this can be done with a precision of the order $N^{-1/2} + N^{1/2}/m_T$ (up to some correction factor) where $m_T$ denotes the mean number of events per point process. In (Liu, 2020), the author complements the analysis started in (Delattre and Fournier, 2016) by investigating the problem of estimating the density of connections p in the same setting with the difference being that now one has access only to the information present in K Hawkes point processes for some large $K \le N$. Under this constraint, (Liu, 2020) proposes an estimator of the parameter p with rate of convergence $K^{-1/2} + N/(K^{1/2}m_T) + N/(Km_T^{1/2})$.

In the present paper, we consider a graphical model which is defined in terms of a system of N interacting $\{0, 1\}$-valued chains $X = \{X_{i,t}, t \in \mathbb{Z}, 1 \le i \le N\}$ and a random matrix $\theta = (\theta_{ij})_{1 \le i,j \le N}$ of i.i.d entries distributed as Ber(p). The event $X_{i,t} = 1$ indicates that the i-th chain sends some signal at time t, and $X_{i,t} = 0$ otherwise. Let $X_t = (X_{1,t}, \dots, X_{N,t})$ be the configuration of the system X at time t. The model is defined as follows. Conditionally on $\theta$, the system X evolves as a stationary Markov chain on the state space $\{0, 1\}^N$ in which the conditional distribution of $X_t$ given that $X_{t-1} = x$ is that of N independent Bernoulli random variables with parameter $p_{\theta,i}(x)$, $i = 1, \dots, N$, where

$$ (1) \quad p_{\theta,i}(x) = \mu + (1 - \lambda) \left( \frac{1}{N} \sum_{j \in \mathcal{P}_+} \theta_{ij} x_j + \frac{1}{N} \sum_{j \in \mathcal{P}_-} \theta_{ij} (1 - x_j) \right), \quad x = (x_1, \dots, x_N). $$

Here, $0 < \lambda < 1$ and $0 \le \mu \le \lambda$ are parameters, and $\mathcal{P}_+$, $\mathcal{P}_-$ are subsets of $[N] = \{1, \dots, N\}$ forming a partition. Neither the parameters $\mu$ and $\lambda$ nor the subsets $\mathcal{P}_+$ and $\mathcal{P}_-$ depend on the random matrix $\theta$.

The function $p_{\theta,i}(x)$ appearing in (1) models the probability of observing a new signal for the i-th chain given that the configuration of the system at the precedent time is $x$. The parameter $\mu$ models the baseline activity of each interacting chain. Note that the probability $p_{\theta,i}(x)$ depends on the past $x$ only through the values $x_j$ for which $\theta_{ij} = 1$, so that the random matrix $\theta$ encodes the interaction among the chains of the system. In light of this, $\theta$ can be considered as the proxy for the graph of conditional dependencies of the model. Note also that each chain is either excitatory or inhibitory in the model; the set $\mathcal{P}_+$ denotes the set of excitatory chains and $\mathcal{P}_-$ the set of inhibitory ones. An excitatory chain increases the probability of future signals to occur whenever it sends any signal, whereas the signals sent
by any inhibitory chain reduce the probability of occurrence of new signals. The increment of this changing is $(1 - \lambda)/N$, implying that the interactions are of mean-field type and that $1 - \lambda$ parametrizes the strength of the interactions.

In what follows, we denote $r_+^N = |\mathcal{P}_+|/N$ and $r_-^N = |\mathcal{P}_-|/N$ the fraction of excitatory and inhibitory components of the model, respectively, and suppose that there are values $0 < r_+, r_- < 1$ satisfying $r_+ + r_- = 1$ such that $|r_+^N - r_+| \vee |r_-^N - r_-| \le KN^{-1}$ for some universal constant $K$. One can easily check that a choice satisfying this assumption with $K=1$ is $|\mathcal{P}_+| = \lceil r_+N \rceil$ and $|\mathcal{P}_-| = N - |\mathcal{P}_+|$.

**Problem formulation.** The goal of this paper is to address the following statistical question. By observing a sample $X_1, \dots, X_T$ of the system $X$, can we estimate the asymptotic density of connections $p$ of the random matrix $\theta$, knowing only the size of the system $N$ and the asymptotic fraction of excitatory and inhibitory components $r_+$ and $r_-$? This means that we do not assume any prior knowledge on $\mu, \lambda, p, \theta, \mathcal{P}_+, \mathcal{P}_-$, i.e., they are all unknown.

Although the estimation problem investigated here is similar to the ones considered in (Delattre and Fournier, 2016; Liu, 2020), our work is different from theirs in at least two substantial ways. First, in their model the interaction between each pair of components can only be excitatory, whereas it can be either excitatory or inhibitory in ours, making the analysis of the estimation problem (especially the study of the random environment) more challenging. Second, a crucial step in the analysis is to study the decay of correlation between $X_{i,t}$ and $X_{j,s}$ (and of products of these), for each fixed realization of the random matrix $\theta$. In our case, this is done through the study of coalescing random walks defining a backward regeneration representation of the system $X$. Interestingly, we show that this backward regeneration representation allows us to perform perfect simulation of the system $X$. That is, one can exactly sample the system $X$ (conditionally on $\theta$) from its unique stationary distribution. These probabilistic results have an interest in its own. In (Delattre and Fournier, 2016) and (Liu, 2020), on the other hand, the analysis of the model for each realization of the graph of conditional dependence follows a rather different approach, fundamentally based on martingale arguments. Moreover, the question of how to perfectly sample the underlying process is not discussed.

As will become apparent in our results, we can not only estimate the parameter $p$ but also the parameters $\mu$ and $\lambda$. In our main result, we show that the parameters $(\mu, \lambda, p)$ can be simultaneously estimated with rate $N^{-1/2} + N^{1/2}/T + (\log(T)/T)^{1/2}$, under some extra conditions on the asymptotic fraction $r_+$ of excitatory components of the model. This convergence rate is consistent with the one found in (Delattre and Fournier, 2016) (up to a $\log(T)$ factor), because the quantity $m_T$ (adapted to our setting) increases linearly in $T$. In Appendix I, we briefly discuss the optimality of our estimation rate by analyzing two related statistical settings in which the observations follow a binomial mixture model. Taken together, the results in Appendix I suggest that our estimation rate may be near optimal. Establishing rigorously such a result in the statistical setting considered in this paper remains an interesting open problem.

We are particularly motivated by the statistical analysis of neuronal networks. Neurons communicate with each other by firing short electrical pulses, often called spikes. The spiking times of a network of neurons depend on its graph of interaction, a combinatorial structure (typically unknown) encoding the type of interaction (excitatory or inhibitory) between each pair of neurons in the network. An important question in Neuroscience is to understand which features of this graph can be inferred from the spiking times of the recorded neurons. Our results suggest that one such feature is the density of connections in the graph, as long as we know the true fraction $r_+$ of excitatory neurons in the network. In many situations, this turns out to be the case. For example, in humans and many other mammalian species the value $r_+$ is known and ranges from 0.7 to 0.8.
Analogously to Delattre and Fournier (2016), our results also hold in the symmetric interaction case, namely when $(\theta_{ij})_{1\le i\le j\le N}$ are i.i.d. and $\theta_{ij} = \theta_{ji}$ for all $1 \le i, j \le N$. Only a few steps in the proofs need to be adapted when working in this new framework, and all these modifications are located in Appendix F. We have gathered the details at the beginning of that section.

Throughout the paper, the parameter $p$ does not scale with the number of components $N$ of the system $X$, ensuring that the random matrix $\theta$ remains typically dense. An interesting open question is whether our results extend to sparse regimes, where $p = c/N$ for some positive constant $c > 0$ or $p \sim \log(N)/N$. In these sparse regimes, the model must be suitably modified, as the normalization factor $N^{-1}$ is no longer appropriate. Furthermore, many of the results concerning the coalescing random walks associated with the backward regeneration representation of the system $X$ (e.g. Proposition 3.5) do not directly apply and require substantial extensions. Such extensions are beyond the scope of this work.

Finally, let us highlight two recent works (Gaitonde, Moitra and Mossel, 2024) and (Chevallier and Ost, 2024), which appeared during the revision of this paper. The work of (Gaitonde, Moitra and Mossel, 2024) deals with the problem of estimating order $k$ graphical models with $N$ sites and bounded degree from dependent trajectories generated by Glauber dynamics. Among their results, the authors show that given $O(N \log N)$ total site updates from the Glauber trajectory, it is possible to correctly output the conditional dependency structure of the underlying graphical model in time $O(N^2 \log N)$, overcoming the known $N^{\Theta(k)}$ computational barrier of the i.i.d. setting. Our framework differs from theirs in two key aspects. First, the dynamics of their model is reversible, whereas ours is not. Second, their analysis assumes that the chosen site to be updated is always observed at each time step, even when the configuration of the system remains unchanged, an assumption that may be difficult to verify in some real datasets. The work of (Chevallier and Ost, 2024) complements the analysis done in the present paper by addressing the problem of estimating the sets $P_+$ and $P_-$ without prior knowledge of the other model parameters. We believe these results could enable us to estimate edge probability $p$ even without knowing the asymptotic fraction of excitatory and inhibitory components $r_+$ and $r_-$ of the system $X$. For sake of simplicity, we do not pursue this extension in the present work.

In the next section we define our model rigorously, introduce some general notation used throughout the paper, state our main results and give some heuristics behind them. At the end of that section, we also provide the organization of the rest of the paper.

## 2. Model definition, notation and main results.

*Model definition.* We consider a system of $N$ interacting chains $X = \{X_{i,t}, t \in \mathbb{Z}, 1 \le i \le N\}$ taking values in $\{0, 1\}$ denoting the presence or the absence of a signal at a given time. This system evolves in a random environment which is given by the realization of a directed Erdös-Rényi random graph via the selection of $N^2$ i.i.d. Bernoulli random variables $\theta_{ij} \sim Ber(p), 1 \le i, j \le N$ with $0 \le p \le 1$. Conditionally on the realization of $\theta = (\theta_{ij})_{1\le i,j\le N}$, the evolution of the system $X$ is that of a stationary Markov chain on the state space $\{0, 1\}^N$ with transition probabilities (which depend on $\theta$) given as follows. Writing $X_t = (X_{1,t}, \dots, X_{N,t})$ and $x, y \in \{0, 1\}^N$, we have, for all $t \in \mathbb{Z}$,

$$ (2) \qquad IP_{\theta}(X_t = y | X_{t-1} = x) = \prod_{i=1}^{N} (p_{\theta,i}(x))^{y_i} (1 - p_{\theta,i}(x))^{(1-y_i)}, $$

where $p_{\theta,i}(x) = IP_{\theta}(X_{i,t} = 1 | X_{t-1} = x)$ is defined in (1). The existence and uniqueness of a stationary version of the Markov chain having transition probabilities as defined in (2) and (1) is granted by Theorem 3.3 presented in Section 3.
**Notation.** Throughout the paper, the letters $t, s$ denote some time values whereas the letters $i, j$ denote some spatial values, that is the index of the corresponding component. For a vector $v$ in $\mathbb{R}^N$, the notation $\bar{v}$ denotes the spatial mean $\bar{v} = N^{-1} \sum_{i=1}^N v_i$. In agreement with the left-hand side of (2), we write $\mathbb{P}_{\theta}$ to denote the probability measure under which the environment $\theta$ is kept fixed and the process $X$ is distributed as the unique stationary version of the Markov chain having transition probabilities as defined in (2) and (1). The expectation taken with respect to $\mathbb{P}_{\theta}$ is denoted by $\mathbb{E}_{\theta}$. The variance and covariance computed from $\mathbb{E}_{\theta}$ are denoted $\mathrm{Var}_{\theta}$ and $\mathrm{Cov}_{\theta}$ respectively. Moreover, we write $\mathbb{P}$ to denote the probability measure under which the random environment $\theta = (\theta_{ij})_{1\le i,j\le N}$ is distributed as a collection of i.i.d. random variables with distribution $\mathrm{Ber}(p)$ and the conditional distribution of the process $X$ given $\theta$ is that of the process $X$ under $\mathbb{P}_{\theta}$, i.e., the following identity holds $\mathbb{P}(X \in \cdot | \theta) = \mathbb{P}_{\theta}(X \in \cdot)$. Finally, we denote $\mathbb{E}$ the expectation taken with respect to the probability measure $\mathbb{P}$, and $\mathrm{Var}$ and $\mathrm{Cov}$ the variance and covariance, respectively, computed from $\mathbb{E}$.

**Main results.** We are given a sample $X_1, \dots, X_T$ of a stationary Markov chain with transition probabilities given by (2) and (1), associated to some realization of the random environment $\theta$ which is not observed. The goal is to estimate the parameter $p$, the asymptotic density of connections in the underlying random environment. We consider an estimator of $p$ which is a function of three other estimators defined below. In what follows, let $\bar{X}_t = N^{-1} \sum_{i=1}^N X_{i,t}$ denote the average number of signals emitted by the system at time $t$; $Z_{i,t} = \sum_{s=1}^t X_{i,s}$ denote the number of signals emitted by the $i$-th component of the system in the discrete interval $\{1, \dots, t\}$ and $\bar{Z}_t = N^{-1} \sum_{i=1}^N Z_{i,t}$ denote the average number of signals emitted by the system in the interval $\{1, \dots, t\}$. For notational convenience, we set $\bar{Z}_0 := 0$. With this notation, our three main estimators are defined as follows:

$$ (3) \qquad \hat{m} = \frac{\bar{Z}_T}{T}, \quad \hat{v} = \frac{(T+1)N}{T^3} \left[ \frac{1}{N} \sum_{i=1}^N (Z_{i,T})^2 - \frac{T}{(T+1)} \left( \bar{Z}_T + (\bar{Z}_T)^2 \right) \right] \quad \text{and} $$

$$ (4) \qquad \hat{w} = 2W_{2\Delta} - W_{\Delta}, \text{ with } W_{\Delta} = \frac{N}{T} \sum_{k=1}^{\lfloor T/\Delta \rfloor} (\bar{Z}_{k\Delta} - \bar{Z}_{(k-1)\Delta} - \Delta\hat{m})^2, $$

where $\Delta \in \{1, \dots, \lfloor T/2 \rfloor\}$ is a tuning parameter. The estimator $\hat{m}$ is the spatio-temporal mean of signals emitted by the system in the observed interval $\{1, \dots, T\}$. As we explain in Section 2.2, the estimator $\hat{v}$ is related to the empirical variance of the random variables $Z_{1,T}, \dots, Z_{N,T}$, i.e., a variance over the spatial components of the system. The estimator $\hat{w}$ is computed from the empirical variance of the random variables $\bar{Z}_{\ell\Delta}, \bar{Z}_{2\ell\Delta} - \bar{Z}_{\ell\Delta}, \dots, \bar{Z}_{T-\ell\Delta}$, i.e., the empirical variance of the mean number of signals emitted by the system over different time intervals. For these reasons, the estimators $\hat{v}$ and $\hat{w}$ are called spatial variance and temporal variance respectively.

Obviously, the three estimators $\hat{m}, \hat{v}, \hat{w}$ depend on $N$ and $T$ (as well as the tuning parameter $\Delta$ for $\hat{w}$). We chose not to specify this dependence in the notation to keep notation as light as possible.

**THEOREM 2.1.** There exists a constant $K > 0$ depending only on $\lambda$ such that for all $\varepsilon \in (0, 1)$, $N \ge 1$, $T \ge 2$ and $1 \le \Delta \le \lfloor T/2 \rfloor$,

$$ (5) \qquad \mathbb{P}(|\hat{m} - m| \ge \varepsilon) \le \frac{K}{\varepsilon} \left( \frac{1}{\sqrt{TN}} + \frac{1}{N} \right), $$
$$ (6) \qquad \mathbb{P}(|\hat{v} - v| \ge \varepsilon) \le \frac{K}{\varepsilon} \left( \frac{\sqrt{N}}{T} + \frac{1}{\sqrt{N}} \right), $$

$$ (7) \qquad \mathbb{P}(|\hat{w} - w| \ge \varepsilon) \le \frac{K}{\varepsilon} \left( \frac{1}{N} + \frac{(1-\lambda)^\Delta}{\Delta} + \sqrt{\frac{\Delta}{T}} \right), $$

where the vector $(m, v, w)$ is given by

$$ (8) \qquad \begin{cases} m = \frac{\mu+(1-\lambda)pr_-}{1-p(1-\lambda)(r_+-r_-)}, \\ v = (1-\lambda)^2 p(1-p)((m-r_-)^2 + r_+r_-), \\ w = m(1-m) \frac{1+4(1-\lambda)^2p^2r_+r_-}{(1-p(1-\lambda)(r_+-r_-))^2}. \end{cases} $$

The proof of Theorem 2.1 is given at the end of Section 4. The constant $K$ depends only on the parameter $\lambda$. Let us mention that this constant diverges when $\lambda \to 0$ (because of the inversion of the linear system, see for instance the bound of Lemma F.3) and $\lambda \to 1$ (because the dependence between the components of the model becomes too weak, see for instance the bounds of Proposition 3.5).

Once the parameters $(m, v, w)$ are estimated from the sample, one wants to deduce some estimators of the parameters $(\mu, \lambda, p)$. Here is how one can proceed.

Let $\Lambda = \{(\mu, \lambda, p) \in (0, 1)^3: 0 < \mu < \lambda\}$ be the open set of admissible parameters. For all $(\mu, \lambda, p) \in \Lambda$, let us define

$$ (9) \qquad D(\lambda, p) = 1 - (1 - \lambda)p(r_+ - r_-) > 0, $$

the denominator that appears in the expressions of $m$ and $w$. Then, for $k = 1, 2, 3$, let $\Psi_k : \Lambda \to \mathbb{R}$ be defined by

$$ (10) \qquad \left\{ \begin{aligned} \Psi_1(\mu, \lambda, p) &= (\mu + (1 - \lambda)pr_-)/D(\lambda, p) \\ \Psi_2(\mu, \lambda, p) &= (1 - \lambda)^2 p(1 - p)[(\Psi_1(\mu, \lambda, p) - r_-)^2 + r_+r_-] \\ \Psi_3(\mu, \lambda, p) &= \Psi_1(\mu, \lambda, p)[1 - \Psi_1(\mu, \lambda, p)][1 + 4(1 - \lambda)^2 p^2 r_+ r_-]/D(\lambda, p)^2. \end{aligned} \right. $$

Finally, let $\Psi : \Lambda \to \mathbb{R}^3$ be defined by the three coordinate functions above so that Equation (8) rewrites as $(m, v, w) = \Psi(\mu, \lambda, p)$. Finally, let us remark that, whatever the value of $r_+$, the image $\Psi(\Lambda)$ is included in $(0, 1) \times (0, \infty)^2$ (see Proposition G.1 stated and proved in Appendix G).

The following proposition gives some information about the inversion of the function $\Psi$, by means of an auxiliary function $\kappa : (0, 1) \times (0, \infty) \to (0, \infty)$ which is defined by

$$ (11) \qquad \kappa(m, w) = (r_+ - r_-)^2 \frac{w}{m(1 - m)}. $$

**PROPOSITION 2.2.** Whatever the value of $r_+$, there exist two explicit functions $\Phi^{(+)}$ and $\Phi^{(-)}$ such that the following results hold.

1. For all $(\mu, \lambda, p) \in \Lambda$, $(\mu, \lambda, p) \in \{\Phi^{(+)} \circ \Psi(\mu, \lambda, p), \Phi^{(-)} \circ \Psi(\mu, \lambda, p)\}$,

2. Moreover, if $r_+ \ge 1/2$ or

$$ \kappa(\Psi_1(\mu, \lambda, p), \Psi_3(\mu, \lambda, p)) \ge 4r_+r_-, \\ \text{then } (\mu, \lambda, p) = \Phi^{(-)} \circ \Psi(\mu, \lambda, p). $$

Proposition 2.2 is deduced as a corollary of Proposition G.3 stated and proved in the Appendix G. Furthermore, the expressions of the functions $\Phi^{(+)}$ and $\Phi^{(-)}$ are also given in Appendix G.
### REMARK 2.3.
Note that the system of equations defining $X$ is not symmetric with respect to the sets $\mathcal{P}_+$ and $\mathcal{P}_-$. One source of asymmetry is the term $N^{-1} \sum_{j \in \mathcal{P}_-} \theta_{ij}$ which appears in (1) and is denoted by $L_i^{N, \bullet-}$ below. This asymmetry partly explains why the inversion of (8) is always feasible when $r_+ \ge 1/2$, but not when $r_+ < 1/2$.

We are now in position to deduce the following corollary of Theorem 2.1.

### COROLLARY 2.4.
Assume that the condition 2 of Proposition 2.2 is satisfied. Then there exists a constant $K > 0$ depending only on $\lambda$ such that for all $\varepsilon \in (0, 1)$, $N \ge 1$, $T \ge 2$ and $1 \le \Delta \le [T/2]$,

$$ (12) \quad \mathbb{P}(\|\Phi^{(-)}(\hat{m}, \hat{v}, \hat{w}) - (\mu, \lambda, p)\|_{\infty} \ge \varepsilon) \le \frac{K}{\varepsilon} \left( \frac{\sqrt{N}}{T} + \frac{1}{\sqrt{N}} + \frac{(1-\lambda)^{\Delta}}{\Delta} + \sqrt{\frac{\Delta}{T}} \right). $$

### PROOF.
By definition of $\Psi$, remark that Theorem 2.1 gives a control of $\|(\hat{m}, \hat{v}, \hat{w}) - \Psi(\mu, \lambda, p)\|_{\infty}$. Hence, the result follows from the facts that $\Phi^{(-)}$ is locally Lipschitz (see Proposition G.2) and $\Phi^{(-)} \circ \Psi(\mu, \lambda, p) = (\mu, \lambda, p)$. $\square$

In practice, a value must be chosen for the tuning parameter $\Delta$. The question of how to choose $\Delta$ with respect to $T$ is discussed in Section 5. In what follows we present some heuristic arguments leading to the convergence stated in Theorem 2.1.

### 2.1. Heuristics for the spatio-temporal mean.
Let us remind that $X$ follows the stationary version of a Markov chain having transition probabilities as defined in (2) and (1), conditionally on the environment $\theta$. Denote $m^N := (m_1^N, \dots, m_N^N) = \mathbb{E}_{\theta}[X_0]$. Notice that $m^N$ depends on the realization of $\theta$, but we omit to specify this dependence to keep concise notation. Throughout this section, the notation $a \approx b$ is used to express the fact that $a$ and $b$ are expected to be close to each other as either $N$ or $T$ are large enough. By ergodicity, we must have as $T \to \infty$,

$$ \frac{1}{T} \sum_{t=1}^{T} X_t \to m^N, \quad \mathbb{P}_{\theta} - \text{a.s.}, $$

so that we expect for $T$ large enough that

$$ (13) \quad \hat{m} = \frac{1}{N} \sum_{i=1}^{N} \frac{1}{T} \sum_{t=1}^{T} X_{i,t} \approx \frac{1}{N} \sum_{i=1}^{N} m_i^N := \overline{m^N}. $$

Hence, to find the limit of $\hat{m}$ one needs to understand the asymptotic behavior of $\overline{m^N}$. To that end, observe that by taking first the expectation in both sides of (1) and then using the stationarity, one can check that, for all $1 \le i \le N$,

$$ (14) \quad m_i^N = \mu + (1 - \lambda) \left( \frac{1}{N} \sum_{j \in \mathcal{P}_+} \theta_{ij} m_j^N + \frac{1}{N} \sum_{j \in \mathcal{P}_-} \theta_{ij} (1 - m_j^N) \right). $$

Denoting $A^N(i, j) = N^{-1}\theta_{ij}$ for all $1 \le i \le N$ and $j \in \mathcal{P}_+$ and $A^N(i, j) = -N^{-1}\theta_{ij}$ for all $1 \le i \le N$ and $j \in \mathcal{P}_-$, we can then write the above system of equations in matrix form as follows:

$$ m^N = \mu 1_N + (1 - \lambda) (A^N m^N - L^{N, \bullet-}), $$
where $1_N$ denotes the $N$-dimensional vector with all entries equal to 1 and $L^{N, \bullet-} = (L_1^{N, \bullet-}, \dots, L_N^{N, \bullet-})$ is given by $L_i^{N, \bullet-} = \sum_{j \in \mathcal{P}_-} A^N(i, j)$. Therefore, we deduce that

$$ (I_N - (1 - \lambda)A^N) m^N = \mu 1_N - (1 - \lambda)L^{N, \bullet-}, $$

where $I_N$ denotes the identity matrix of size $N$. Let us denote $Q^N = (I_N - (1 - \lambda)A^N)^{-1}$. The random matrix $Q^N$ is well-defined whenever $\lambda > 0$ (see Appendix F). Granted that $Q^N$ is well-defined, we then have that

$$ (15) \qquad m^N = \mu Q^N 1_N - (1 - \lambda) Q^N L^{N, \bullet-}. $$

Then, remark that $L^{N, \bullet-} = A^N 1_{\mathcal{P}_-}$, where $1_{\mathcal{P}_-}$ is the $N$-dimensional vector with value 1 for the coordinates in $\mathcal{P}_-$ and value 0 otherwise. Using the series expansion $Q^N = \sum_{k=0}^{\infty} (1 - \lambda)^k (A^N)^k$, one can deduce that $Q^N A^N = (Q^N - I_N)/(1 - \lambda)$. Hence, $(1 - \lambda)Q^N L^{N, \bullet-} = \ell^{N, \bullet-} - 1_{\mathcal{P}_-}$ with $\ell^{N, \bullet-} = Q^N 1_{\mathcal{P}_-}$. Similarly, we denote $\ell^N = Q^N 1_N$, so that Equation (15) rewrites as

$$ m^N = \mu \ell^N + 1_{\mathcal{P}_-} - \ell^{N, \bullet-}. $$

Now, let us give some heuristics on the asymptotic behavior of $\ell^N$ and $\ell^{N, \bullet-}$. To that end, denote $L^N = A^N 1_N$ and observe that by the definition of $\ell^N$, using once more the series expansion of $Q^N$,

$$ \ell^N = 1_N + (1 - \lambda)L^N + \sum_{k=2}^{\infty} (1 - \lambda)^k (A^N)^k 1_N. $$

Now, by law of large numbers, we expect that, whenever $N$ is large enough, $L_i^N = \sum_{j=1}^N A^N(i, j) \approx p(r_+ - r_-)$ for each $1 \le i \le N$, and by induction, we expect that

$$ \left( (A^N)^k 1_N \right)_i \approx [p(r_+ - r_-)]^{k-1} L_i^N, $$

for each $1 \le i \le N$. Hence, for $N$ sufficiently large

$$ \ell^N \approx 1_N + \frac{1}{D(\lambda, p)} (1 - \lambda) L^N, $$

where $D(\lambda, p)$ is defined in (9). Similarly, by definition of $\ell^{N, \bullet-}$, we have

$$ \ell^{N, \bullet-} = 1_{\mathcal{P}_-} + (1 - \lambda)L^{N, \bullet-} + \sum_{k=2}^{\infty} (1 - \lambda)^k (A^N)^k 1_{\mathcal{P}_-}. $$

Using the law of large numbers once more, we expect that $L^{N, \bullet-} \approx -pr_- 1_N$ whenever $N$ is large enough, so that by the previous induction, we expect that

$$ (A^N)^k 1_{\mathcal{P}_-} \approx -pr_- [p(r_+ - r_-)]^{k-2} L^N. $$

Hence,

$$ \ell^{N, \bullet-} \approx 1_{\mathcal{P}_-} + (1 - \lambda)L^{N, \bullet-} - \frac{(1 - \lambda)^2 pr_-}{D(\lambda, p)} L^N. $$

All in all, we expect that

$$ (16) \qquad m^N \approx \mu \left( 1_N + \frac{1 - \lambda}{D(\lambda, p)} L^N \right) - (1 - \lambda)L^{N, \bullet-} + \frac{(1 - \lambda)^2 pr_-}{D(\lambda, p)} L^N. $$
Now, remind that $L^N \approx p(r_+ - r_-)1_N$ and $L^{N, \bullet -} \approx -pr_-1_N$, so that

$$ (17) \qquad \overline{m^N} \approx \mu \left( 1 + \frac{(1-\lambda)p(r_+ - r_-)}{D(\lambda, p)} \right) + (1-\lambda)pr_- + \frac{(1-\lambda)^2 p^2 (r_+ - r_-)r_-}{D(\lambda, p)} = m, $$

where m is defined in (8). Combining (13) and (17), we expect that for T and N large enough,

$$ \hat{m} \approx \frac{\mu + (1 - \lambda)pr_{-}}{D(\lambda, p)} = m. $$

This heuristics is made precise in Equations (32) and (35).

## 2.2. Heuristics for the spatial variance.
First note that we can always write that

$$ \mathrm{Var}(Z_{i,T}) - \mathrm{IE}[\mathrm{Var}_{\theta}(Z_{i,T})] = \mathrm{Var}(\mathrm{IE}_{\theta}[Z_{i,T}]). $$

The three terms above are considered from left to right.

First, since the spatial covariance of the model vanishes sufficiently fast (see Lemma 3.7 for a precise statement), the term $\mathrm{Var}(Z_{i,T})$ can be estimated from the data via an empirical variance, that is,

$$ (18) \qquad \mathrm{Var}(Z_{i,T}) \approx \frac{1}{N} \sum_{i=1}^{N} (Z_{i,T} - \bar{Z}_{T})^2. $$

Second, the term $\mathrm{IE}[\mathrm{Var}_{\theta}(Z_{i,T})]$ can also be estimated from the data. Indeed, we expect that the temporal covariance of the model vanishes sufficiently fast (see Lemma 3.7 for a precise statement) so that, for T large enough,

$$ (19) \qquad \mathrm{IE}[\mathrm{Var}_{\theta}(Z_{i,T})] \approx T \mathrm{IE}[\mathrm{Var}_{\theta}(X_{i,0})]. $$

Given the random matrix $\theta$, we know that $X_{i,0}$ is a Bernoulli variable with parameter $m_i^N$ by stationarity of the process, so that

$$ \mathrm{Var}_{\theta}(X_{i,0}) = m_i^N - (m_i^N)^2. $$

Yet, we expect that $T^{-1}Z_{i,T} \approx m_i^N$ which implies that we can estimate $\mathrm{IE}[\mathrm{Var}_{\theta}(X_{i,0})]$ and get

$$ (20) \qquad \frac{1}{N} \sum_{i=1}^{N} \frac{Z_{i,T}}{T} - \left( \frac{Z_{i,T}}{T} \right)^2 = \frac{\bar{Z}_{T}}{T} - \frac{1}{N} \sum_{i=1}^{N} \left( \frac{Z_{i,T}}{T} \right)^2 \approx \mathrm{IE}[\mathrm{Var}_{\theta}(X_{i,0})]. $$

Finally, we may express the limit of the last term. Since $\mathrm{IE}_{\theta}[Z_{i,T}] = Tm_i^N$, we can start from the heuristics (16) and look at its variance. It is easy to check that $\mathrm{Var}(L_i^N) = N^{-1}p(1-p)$ and $\mathrm{Var}(L_i^{N, \bullet -}) = \mathrm{Cov}(L_i^N, L_i^{N, \bullet -}) = N^{-1}p(1-p)r_-$. Hence, we expect that

$$ \begin{aligned} \mathrm{Var}(m_i^N) &\approx N^{-1}p(1-p) \left[ \left( \frac{\mu(1-\lambda) + (1-\lambda)^2pr_-}{D(\lambda, p)} \right)^2 - 2 \left( \frac{\mu(1-\lambda) + (1-\lambda)^2pr_-}{D(\lambda, p)} \right) (1-\lambda)r_- + (1-\lambda)^2r_- \right] \\ &= N^{-1}(1-\lambda)^2p(1-p)(m^2 - 2mr_- + r_-) = N^{-1}v, \end{aligned} $$

where m and v are defined in (8) and $D(\lambda, p)$ is defined in (9). Therefore,

$$ (21) \qquad \frac{N}{T^2} \mathrm{Var}(\mathrm{IE}_{\theta}[Z_{i,T}]) \approx v. $$
Combining Equations (18)-(21), we have

$$ \hat{v} = \frac{1}{T^2} \sum_{i=1}^{N} (Z_{i,T} - \bar{Z}_T)^2 - \frac{N}{T^2} \bar{Z}_T + \frac{1}{T^3} \sum_{i=1}^{N} (Z_{i,T})^2 \approx v. $$

This heuristics is made precise in Equations (33) and (36).

2.3. Heuristics for the temporal variance. Let us remark that

$$ \frac{\Delta}{N} W_{\Delta} = \frac{\Delta}{T} \sum_{k=1}^{T/\Delta} [(\bar{Z}_{k\Delta} - \bar{Z}_{(k-1)\Delta}) - \Delta T^{-1} \bar{Z}_T]^2. $$

Because we expect that the temporal covariance of the model vanishes sufficiently fast, we should have (supposing $\Delta \ll T$) that for $T$ large,

$$ (22) \qquad \frac{\Delta}{N} W_{\Delta} \approx \mathrm{Var}_{\theta}[\bar{Z}_{\Delta}] = \mathrm{Var}_{\theta}[\bar{U}_{\Delta}], $$

where $U_t = (U_{1,t}, \dots, U_{N,t})$ for each $t \ge 1$, with

$$ (23) \qquad U_{i,t} = Z_{i,t} - \mathrm{IE}_{\theta}[Z_{i,t}] \quad \text{and} \quad \bar{U}_t = N^{-1} \sum_{i=1}^{N} U_{i,t}. $$

Next, (1) implies that for all $1 \le i \le N$,

$$ U_{i,t} = M_{i,t} + (1 - \lambda) \sum_{j=1}^{N} A^{N}(i,j) \sum_{s=1}^{t} (X_{j,s-1} - \mathrm{IE}_{\theta}[X_{j,s-1}]), $$

where $M_{i,t} = \sum_{s=1}^{t} (X_{i,s} - p_{\theta,i}(X_{s-1}))$ is a martingale. From the above identity, one can deduce that

$$ (24) \qquad U_t = M_t + (1 - \lambda) A^N U_t + R_t^N, $$

where $M_t = (M_{1,t}, \dots, M_{N,t})$ and $R_t^N = (1 - \lambda) A^N (X_0 - X_t)$. From our analysis, we have that the term $R_t^N$ is negligible (see Lemma E.2), so that it follows from (24) that $U_t \approx M_t + (1 - \lambda) A^N U_t$, or equivalently,

$$ U_t \approx Q^N M_t. $$

Hence, introducing the vector $(Q^N)^{\mathrm{T}} 1_N = c^N = (c_1^N, \dots, c_N^N)$ where $(Q^N)^{\mathrm{T}}$ denotes the transpose of $Q^N$, we deduce that

$$ \bar{U}_t \approx \overline{Q^N M_t} = \frac{1}{N} \sum_{i=1}^{N} \sum_{j=1}^{N} Q^N(i,j) M_{j,t} = \frac{1}{N} \sum_{j=1}^{N} c_j^N M_{j,t}, $$

which, in turn, implies that

$$ \mathrm{Var}_{\theta}[\bar{U}_t] \approx \mathrm{Var}_{\theta}[\overline{Q^N M_t}] = \frac{1}{N^2} \sum_{i=1}^{N} \sum_{j=1}^{N} c_i^N c_j^N \mathrm{IE}_{\theta}[M_{i,t} M_{j,t}]. $$

Since the martingales $(M_{i,t})_t$ and $(M_{j,t})_t$ are orthogonal for any $i \ne j$, i.e. $\mathrm{IE}_{\theta}[M_{i,t} M_{j,t}] = 0$, one obtains that

$$ (25) \qquad \mathrm{Var}_{\theta}[\bar{U}_t] \approx \frac{1}{N^2} \sum_{i=1}^{N} (c_i^N)^2 \mathrm{IE}_{\theta}[(M_{i,t})^2]. $$
Then, by means of the quadratic variation and the stationarity of the process, we have

$$ (26) \qquad \mathbb{E}_{\theta} \left[ (M_{i,t})^2 \right] = \mathbb{E}_{\theta} \left[ \sum_{s=1}^{t} (X_{i,s} - p_{\theta,i}(X_{s-1}))^2 \right] = t \mathbb{E}_{\theta} \left[ (X_{i,1} - p_{\theta,i}(X_0))^2 \right]. $$

Observing that, conditionally on $\theta$, $X_{i,1}$ is distributed as a Bernoulli random variable with parameter $m_i^N$, and using the fact that $p_{\theta,i}(X_0) \approx m_i^N$, since the spatial covariances are weak, we can deduce (see Lemma E.3) that

$$ (27) \qquad \mathbb{E}_{\theta} \left[ (X_{i,1} - p_{\theta,i}(X_0))^2 \right] \approx m_i^N - (m_i^N)^2. $$

Combining (22) and (25)-(27), we expect that

$$ W_{\Delta} \approx \frac{1}{N} \sum_{i=1}^{N} (c_{i}^{N})^{2} (m_{i}^{N} - (m_{i}^{N})^{2}). $$

Now, as $N$ is large, we have $m_i^N - (m_i^N)^2 \approx m(1-m)$ for all $1 \le i \le N$. Finally, to get the asymptotics of $N^{-1} \sum_{i=1}^N (c_i^N)^2$, one can follow the same lines as in Section 2.1 to check that

$$ c^N \approx 1_N + \frac{(1-\lambda)}{D(\lambda, p)} C^N, $$

with $C^N = (A^N)^{\top} 1_N$ and $D(\lambda, p)$ as defined in (9). Then, by observing that $C_k^N \approx p$ for $k \in \mathcal{P}_+$ and $C_k^N \approx -p$ for $k \in \mathcal{P}_-$, one can check that $(c_k^N)^2 \approx (1+a)^2$ for $k \in \mathcal{P}_+$ and $(c_k^N)^2 \approx (1-a)^2$ for $k \in \mathcal{P}_-$, with

$$ a = \frac{(1 - \lambda)p}{D(\lambda, p)}. $$

Hence,

$$ \frac{1}{N} \sum_{k=1}^{N} (c_k^N)^2 \approx 1 + a^2 + 2a(r_+ - r_-) \\ = \frac{1 + (1 - \lambda)^2 p^2 (1 - (r_+ - r_-)^2)}{(D(\lambda, p))^2} = \frac{1 + 4(1 - \lambda)^2 p^2 r_+ r_-}{(D(\lambda, p))^2}. $$

All in all, we find that

$$ W_{\Delta} \approx m(1-m) \frac{1+4(1-\lambda)^2 p^2 r_+ r_-}{(D(\lambda, p))^2} = w. $$

This heuristics is made precise in Equations (34) and (37).

By relying on the above heuristics now replacing $2\Delta$ by $\Delta$, we also deduce that $W_{2\Delta}$ should converge to $w$ so that

$$ \hat{w} = 2W_{2\Delta} - W_{\Delta} \approx w. $$

The reason of estimating $w$ using $\hat{w}$ instead of $W_{\Delta}$ (or $W_{2\Delta}$) alone is to improve the convergence rate (see the beginning of Appendix E.2).
*Organization of the rest of the paper.* The rest of this paper is devoted to the proof of Theorem 2.1. Its proof is an immediate consequence of Propositions 4.1 and 4.2 and is presented at the end of Section 4. The proof of Proposition 4.1 is given in Appendix E. It relies on a precise study of the decay of correlation of $X_{i,t}$ and $X_{j,s}$ (and of products of these) conditionally on a fixed realization of the environment $\theta$. This is done via a backward regeneration scheme presented in Section 3. The proof of Proposition 4.2 is given in Appendix F. It relies on a study of the random matrices $A^N$ and $Q^N$ and their associated row and column sums. In Section 5, we illustrate the performance of our estimators through simulations. Finally, many other auxiliary results used throughout the paper are proved in appendices.

**3. Backward regeneration scheme.** Recall formula (1). First of all, since $\mu \le \lambda$, let us introduce $\beta := \mu/\lambda$ such that $0 \le \beta \le 1$. Then we may interpret formula (1) in the following way. At any given time $t$, the $i$-th component first decides to update independently of anything else with probability $\lambda$. If it does so, then it decides to send a signal ($X_{i,t} = 1$) with probability $\beta$, else it does not send a signal. Moreover, if the component does not update independently of anything else, then it chooses one of its $N$ neighbors (including $i$ itself) randomly according to the uniform distribution. Suppose it has chosen $j$ as its neighbor. Then there are three possibilities:

* if $\theta_{ij} = 0$, then $X_{i,t} = 0$,

* if $\theta_{ij} = 1$ and $j \in \mathcal{P}_+$, then $X_{i,t}$ copies the value of $X_{j,t-1}$,

* if $\theta_{ij} = 1$ and $j \in \mathcal{P}_-$, then $X_{i,t}$ copies $1 - X_{j,t-1}$.

This interpretation is formalized in Equation (29) below, where the process is constructed via a backward regeneration scheme.

Models like this are called *imitation models* and have been studied in the one-dimensional frame in (De Santis and Piccioni, 2015). Since $\lambda > 0$, at each time step the probability of making an update which is independent of anything else is strictly positive, and this implies the existence of a unique invariant measure of the system as well as the possibility of perfectly sampling from it, which is formalized in the following section.

### 3.1. Backward regeneration representation.

*Notation.* In what follows, we denote by $z = (i, t)$ some space-time coordinate in $\mathcal{Z} = \{1, \dots, N\} \times \mathbb{Z}$ and we write $X_z$ instead of $X_{i,t}$ for all $z = (i, t) \in \mathcal{Z}$.

To construct the Markov chain $X = \{X_{i,t}, t \in \mathbb{Z}, 1 \le i \le N\}$ via a backward regeneration representation, we introduce the following random variables. To each space-time coordinate $z \in \mathcal{Z}$, we associate a couple of independent random variables $(J_z, \xi_z)$ taking values in $\{0, 1, \dots, N\} \times \{0, 1\}$ such that $IP(J_z = 0) = \lambda$, $IP(J_z = k) = (1 - \lambda) \frac{1}{N}$ for $1 \le k \le N$, and $IP(\xi_z = 1) = 1 - IP(\xi_z = 0) = \beta$. Moreover, we suppose that the couples $((J_z, \xi_z), z \in \mathcal{Z})$ are independent. By convention, let us define $J_{(0,t)} = 0$ for all $t \in \mathbb{Z}$.

Then, for any $z = (i, t) \in \mathcal{Z}$, let us define a backward random walk $I^z = (I_s^z)_{s \in \mathbb{Z}}$ taking values in the state space $\{0, 1, \dots, N, \infty\}$, where $I_s^z = \infty$ for all $s > t$, $I_t^z = i$ and, as $s$ decreases, $I^z$ follows the space coordinates given by the $J$ variables, that is

$$ (28) \qquad I_{s-1}^z = J_{(I_s^z, s)}, \quad \text{for all } s \le t. $$

Because of the convention $J_{(0,t)} = 0$, the state 0 is a cemetery for the process $I^z$. Let us denote $\tau_z^R$ the last time that the random walk $I^z$ is not in the cemetery state 0:

$$ \tau_z^R = \inf\{s \in \mathbb{Z} : I_s^z > 0\}. $$

Reaching the cemetery happens after a finite time, almost surely for all the sites.
**PROPOSITION 3.1.** For any $\theta$, it holds that $\mathbb{P}_{\theta}(\forall z \in \mathcal{Z}, \tau_z^R > -\infty) = 1$.

**PROOF.** Since $\mathcal{Z}$ is countable, it suffices to check that $\mathbb{P}_{\theta}(\tau_z^R > -\infty) = 1$ for all $z \in \mathcal{Z}$.
For each $k \ge 1$,

$$ \mathbb{P}(\tau_z^R \le t - k) = \mathbb{P}(I_{t-k}^z \neq 0, \dots, I_{t-1}^z \neq 0) \le (1 - \lambda)^k, $$

where we used the independence between the random variables $J_{\bar{z}}$ and the fact that $\mathbb{P}(J_{\bar{z}} \neq 0) = 1 - \lambda$. It follows that $\mathbb{P}(\tau_z^R = -\infty) = \lim_{k \to \infty} \mathbb{P}(\tau_z^R \le t - k) = 0$. $\square$

The time $\tau_z^R$ is called a regeneration time for state $z$. Let us then define the set of regenerating sites by

$$ \mathcal{R} = \{z \in \mathcal{Z} : J_z = 0\}. $$

**REMARK 3.2.** The sites $z = (i, t)$ such that $\theta_{iJ_z} = 0$ and $J_z \ge 1$ could also be considered as regenerating sites. Since we are interested in results valid for any realization of $\theta$, we did not include those sites in the definition of $\mathcal{R}$.

Observe that since $I_{\tau_z^R - 1}^z = J_{(\tau_z^R, \tau_z^R)} = 0$, we have that $(I_{\tau_z^R, \tau_z^R}^z) \in \mathcal{R}$ is a regenerating site.

In summary, the random walk $I^z$ equals $+\infty$ for all times $u > t$, that is, before it “starts to live”. Then, the random walk lives in $\{1, \dots, N\}$ until it reaches 0 and thus regenerates. It then remains in state 0 forever.

For all $z = (i, t) \in \mathcal{Z}$, we define $X_z$ as follows:

$$ (29) \qquad X_z = \xi_z \mathbf{1}_{z \in \mathcal{R}} + \theta_{iJ_z} (X_{J_z, t-1} \mathbf{1}_{J_z \in \mathcal{P}_+} + (1 - X_{J_z, t-1}) \mathbf{1}_{J_z \in \mathcal{P}_-}) \mathbf{1}_{z \notin \mathcal{R}}. $$

Note that, if $z \in \mathcal{R}$, then $X_z$ is sampled according to the Bernoulli random variable $\xi_z$. When $z \notin \mathcal{R}$, the random variables $X_z$ are defined recursively. On $\Omega_0 = \{\forall z \in \mathcal{Z}, \tau_z^R > -\infty\}$, this recursion ends in finite time for every site and so the process $X$ constructed via (29) is well-defined. Proposition 3.1 ensures that $\mathbb{P}_{\theta}(\Omega_0) = 1$, implying that the process $X$ is well-defined almost surely. Since the process $X$ is defined in (29) via the i.i.d. couples $(J_z, \xi_z)$, one can prove that it is a stationary version of the studied Markov chain.

**THEOREM 3.3.** Let $X = (X_z)_{z \in \mathcal{Z}}$ be the process defined through (29). Then, conditionally on $\theta$, $X$ is a stationary Markov chain with transition probability given by (2) and (1).

The proof of Theorem 3.3 is given in Appendix A.

**REMARK 3.4.** The representation (29) is an example of a more general construction that allows to perfectly sample from a process using a *clan of ancestors method*. This method works under certain structural conditions on the transition probabilities of the process; that is, supposing that it is possible to decompose the transition probabilities into a mixture of more elementary transitions. These ideas go back at least to Ferrari (1990), followed by Ferrari et al. (2000) and Comets, Fernández and Ferrari (2002). We also refer to Fernández, Ferrari and Galves (2001) for a comprehensive introduction to this approach.

The regeneration representation presented in (29) above induces a natural (random) partitioning of the state space $\mathcal{Z}$ via the equivalence relation of coalescence defined by

$$ z_1 \rightsquigarrow z_2 \text{ if and only if } \exists s \in \mathbb{Z}, I_s^{z_1} = I_s^{z_2} \notin \{0, \infty\}. $$

The next section is devoted to a study of this equivalence relation.
### 3.2. Coalescence of the backward random walks.
For any $z_1, z_2 \in \mathcal{Z}$, let us define the random time

$$ (30) \qquad \tau_{z_1, z_2}^c = \sup \{s \in \mathbb{Z} : I_s^{z_1} = I_s^{z_2} \notin \{0, \infty\}\}, $$

namely the time at which $z_1$ and $z_2$ coalesce (with the convention that $\tau_{z_1, z_2}^c = -\infty$ if they do not coalesce, denoted by $z_1 \rightsquigarrow z_2$). The coalescence property is related to the dependence between the coordinates of the process $X$. Loosely speaking, if $z_1$ does not coalesce with $z_2$ then $X_{z_1}$ behaves as a copy of itself which is independent of $X_{z_2}$. This idea is formalized via a coupling argument in Appendix B.

In the next result, we provide upper bounds for the probability of events involving the coalescence of two or more backward random walks.

**PROPOSITION 3.5.** Let $z_k = (i_k, t_k) \in \mathcal{Z}$, for $k = 1, \dots, 4$, be four different sites. There exists a constant $K$ only depending on $\lambda$ such that the following inequalities hold.

(i) $\mathbb{P}_{\theta}(\{z_1 \rightsquigarrow z_2\}) \le \frac{(1-\lambda)^{|t_1-t_2|\vee 1}}{1-(1-\lambda)^2} N^{-1}$.

(ii) $\mathbb{P}_{\theta}(\{z_1 \rightsquigarrow z_2 \rightsquigarrow z_3\}) \le K(1-\lambda)^{\langle t-t \rangle} N^{-2}$, where $\bar{t} = t_1 \vee t_2 \vee t_3$ and $\underline{t} = t_1 \wedge t_2 \wedge t_3$.

(iii) $\mathbb{P}_{\theta}(\{z_1 \rightsquigarrow z_2\} \cap \{z_3 \rightsquigarrow z_4\}) \le K \left( (1-\lambda)^{|t_1-t_2|+|t_3-t_4|} N^{-2} + (1-\lambda)^{\langle t-t \rangle} N^{-3} \right)$, where $\bar{t} = t_1 \vee t_2 \vee t_3 \vee t_4$ and $\underline{t} = t_1 \wedge t_2 \wedge t_3 \wedge t_4$.

(iv) If $t_1 \wedge t_2 \ge t_3 \vee t_4$,

$$ \mathbb{P}_{\theta}(\{z_1 \rightsquigarrow z_3 \rightsquigarrow z_4 \not\rightsquigarrow z_2\} \cap \{\tau_{z_2}^R < t_3 \vee t_4\}) \le K(1-\lambda)^{t_1+t_2-t_3-t_4} N^{-2}. $$

The proof of Proposition 3.5 is given in Appendix C.

**REMARK 3.6.** Here are two simple remarks to intuitively understand the upper bounds of Proposition 3.5.

1. The exponent of the factor $(1 - \lambda)$ is the minimal waiting time before coalescence. For instance, in Item (ii), the first-born site (with birth time $\bar{t}$, remind that time runs backwards) must wait (hence survive) at least $\bar{t} - \underline{t}$ time steps before coalescing with the last-born site.

In that respect, Item (iv) is a bit different because of the event $\{\tau_{z_2}^R < t_3 \vee t_4\}$ which enforces the site $z_2$ to survive even if it does not coalesce.

2. Each coalescence event adds a factor $N^{-1}$. For instance, in Items (ii)-(iv), at least two coalescences are needed.

For all $z \in \mathcal{Z}$, let us denote $Y_z = X_z - \mathbb{E}_{\theta}(X_z)$ the centered version of $X_z$ and observe that $|Y_z| \le 1$ almost surely. Since (the absence of) coalescence is related with independence, it is natural to expect that the rates of decay of coalescence probabilities in Proposition 3.5 can be transferred to rates of decay of covariances for the process. Indeed, using Proposition 3.5 and the coupling introduced in Appendix B, we obtain several upper bounds on the covariance of the variables $Y_z$'s (Lemma 3.7) and the products of these variables (Lemma 3.8). This last result is the main result of this section.

**LEMMA 3.7.** Let $z_1 = (i_1, t_1)$ and $z_2 = (i_2, t_2)$ two points in $\mathcal{Z}$. There exists a constant $K$ only depending on $\lambda$ such that: if $z_1 \ne z_2$,

$$ |\mathrm{Cov}_{\theta}[X_{z_1}, X_{z_2}]| = |\mathrm{Cov}_{\theta}[Y_{z_1}, Y_{z_2}]| \le K(1-\lambda)^{|t_1-t_2|\vee 1} N^{-1}. $$

Otherwise, the quantity above is obviously bounded by 1.
PROOF. The proof is the same as the one for item 3 of Lemma 3.8 below. □

LEMMA 3.8. Let $z_k = (i_k, t_k) \in \mathcal{Z}$, $k = 1, \dots, 4$, and denote $B = \text{Cov}_{\theta}[Y_{z_1} Y_{z_2}, Y_{z_3} Y_{z_4}]$ and $E = \{z_k : k = 1, \dots, 4\}$. There exists a constant $K$ only depending on $\lambda$ such that:

1. If $\#E = 1$, then $|B| \le 1$.

2. If $\#E = 2$, $z_1 \neq z_2$ and $z_3 \neq z_4$, then $|B| \le 1$.

3. If $\#E = 2$, $z_1 = z_2$ and $z_3 = z_4$, then

$$|B| \le K(1-\lambda)^{|t_1-t_3|} N^{-1}.$$

4. If $\#E = 3$ and $z_1 = z_2$ or $z_3 = z_4$, then

$$|B| \le K(1-\lambda)^{|s_3-s_2|+|s_2-s_1|} N^{-2},$$

where $s_1 \le s_2 \le s_3$ is an ordering of $\{t_1, t_2, t_3, t_4\}$.

5. If $\#E = 3$, $z_1 \neq z_2$ and $z_3 \neq z_4$, for instance assume that $z_1 = z_3$ (and so $z_2 \neq z_4$), then

$$|B| \le K(1-\lambda)^{|t_2-t_4|} N^{-1}.$$

6. If $\#E = 4$ and $t_1 \wedge t_2 \ge t_3 \vee t_4$, then

$$|B| \le K(1-\lambda)^{t_1 \vee t_2 - t_3 \wedge t_4} N^{-2}.$$

7. If $\#E = 4$, then

$$|B| \le K(1-\lambda)^{|s_2-s_1|+|s_3-s_4|} N^{-2},$$

where $s_1 \le \dots \le s_4$ is an ordering of $\{t_1, \dots, t_4\}$.

REMARK 3.9. Observe that we cannot have $z_1 = z_2$ and $z_3 = z_4$ simultaneously in Item 4 of Lemma 3.8. In particular, in the case $z_1 = z_2$ (resp. $z_3 = z_4$), the set of time indices $\{t_1, t_2, t_3, t_4\}$ reduces to the set $\{t_1, t_3, t_4\}$ (resp. $\{t_1, t_2, t_3\}$).

The proof of Lemma 3.8 is given in Appendix D.

**4. Key steps and proof of Theorem 2.1.** In order to quantify the convergence stated in Theorem 2.1, we chose first to study the convergence of our estimators as $T \to \infty$ for a fixed environment and then to study the convergence of these temporal limits with respect to the random environment as $N \to \infty$.

Hereafter, $1_N$ denotes the $N$ dimensional vector full of ones and, for any vector $v \in \mathbb{R}^N$, $\bar{v} = N^{-1} \sum_{i=1}^N v_i$ denotes the arithmetic mean of the coordinates of $v$. Remind our three estimators $\hat{m}$, $\hat{v}$, $\hat{w}$ defined in (3)-(4). Their respective limits when $T \to \infty$ while $\theta$ is fixed, as proved below, are:

$$ (31) \qquad \begin{cases} m_{\infty}^N := \overline{m^N} = \mu \overline{\ell^N} - (1-\lambda) \overline{Q^N L^{N,\bullet}}, \\ v_{\infty}^N := \|m^N - m_{\infty}^N 1_N\|_2^2, \\ w_{\infty}^N := \frac{1}{N} \sum_{i=1}^N (c_i^N)^2 (m_i^N - (m_i^N)^2), \end{cases} $$

where $m^N$ is defined in Section 2.1 and the second equality of the first line comes from Equation (15). The rates of these temporal convergences are given below.
**PROPOSITION 4.1.** There exists a constant $K$ depending only on $\lambda$ such that for all $T \ge 1$, $N \ge 1$,

$$ (32) \qquad \mathbb{E}_{\theta} [ |\hat{m} - m_{\infty}^{N}|^2 ] = \mathrm{Var}_{\theta}(\hat{m}) \le K(TN)^{-1}, $$

$$ (33) \qquad \mathbb{E}_{\theta} [ |\hat{v} - v_{\infty}^{N}| ] \le K \left( 1 + \sqrt{v_{\infty}^{N}} \right) \left( \frac{N}{T^2} + \frac{1}{\sqrt{T}} + \frac{\sqrt{N}}{T} \right). $$

$$ (34) \qquad \mathbb{E}_{\theta} [ |\hat{w} - w_{\infty}^{N}| ] \le K \left( \frac{\Delta}{T} + \frac{1}{N} + \frac{(1-\lambda)^{\Delta}}{\Delta} \right), $$

The proof of Proposition 4.1 is given in Appendix E.

The last ingredient needed to prove Theorem 2.1 are the rates of convergence of the temporal limits $m_{\infty}^N, v_{\infty}^N, w_{\infty}^N$ as $N \to \infty$.

**PROPOSITION 4.2.** Assume that $1 \ge \lambda > 0$. There exists a constant $K > 0$ which depends on $\lambda$ such that for all $N \ge 1$, it holds that

$$ (35) \qquad \mathbb{E} [ |m_{\infty}^{N} - m|^{2} ] \le \frac{K}{N^{2}}, $$

$$ (36) \qquad \mathbb{E} [ |v_{\infty}^{N} - v| ] \le \frac{K}{\sqrt{N}}, $$

$$ (37) \qquad \mathbb{E} [ |w_{\infty}^{N} - w|^{2} ] \le \frac{K}{N^{2}}, $$

where $m, v, w$ are defined in (8)

The proof of Proposition 4.2 is given in Appendix F.

Finally, the proof of Theorem 2.1 is merely the combination of Propositions 4.1 and 4.2.

**PROOF OF THEOREM 2.1.** Let us prove Inequality (6). The proof of the two other inequalities are similar and even simpler (the only subtlety is to remark that $\Delta/T < \sqrt{\Delta/T}$ for Inequality (7)).

First, observe that (6) holds automatically when $N^{1/2}/T > 1$. Suppose now that $N^{1/2}/T^2 \le 1$. By the triangle inequality, we have

$$ \mathbb{E} [ |\hat{v} - v| ] \le \mathbb{E} [ |\hat{v} - v_{\infty}^{N}| ] + \mathbb{E} [ |v_{\infty}^{N} - v| ]. $$

By the sub-additivity of $x \mapsto \sqrt{x}$,

$$ (38) \qquad \mathbb{E} \left[ \sqrt{v_{\infty}^{N}} \right] \le \left( \mathbb{E} \left[ \left\| m^{N} - \overline{m}^{N} 1_{N} \right\|_{2}^{2} \right] \right)^{1/2} \le \left( \mathbb{E} \left[ \left\| m^{N} - \overline{m}^{N} 1_{N} \right\|_{2}^{2} - v \right] \right)^{1/2} + \sqrt{v} \le K. $$

Thus, integrating (33) against $\mathbb{P}$ and using the upperbound above, we obtain

$$ \mathbb{E} [ |\hat{v} - v_{\infty}^{N}| ] \le K \left( \frac{N}{T^2} + \frac{1}{\sqrt{T}} + \frac{\sqrt{N}}{T} \right). $$
In turn, thanks to (36), we have

$$ \mathbb{E}[|\hat{v} - v|] \le K \left( \frac{N}{T^2} + \frac{1}{\sqrt{T}} + \frac{\sqrt{N}}{T} + \frac{1}{\sqrt{N}} \right). $$

Combining this with Markov inequality and then using that $N/T^2 \le N^{1/2}/T$ and $T^{-1/2} = N^{1/4}T^{-1/2}N^{-1/4} \le N^{1/2}T^{-1} + N^{-1/2}$, we obtain that for any $\epsilon \in (0, 1)$,

$$ \mathbb{P}(|\hat{v} - v| \ge \epsilon) \le \frac{K}{\epsilon} \left( \frac{N^{1/2}}{T} + N^{-1/2} \right), $$

proving (6).

<figure><img src="image_3.png" /></figure>

**5. Simulation study.** Before showing some numerical results, the practical implementation of the method (especially the inversion of the function $\Psi$ defined by (10)) and the choice of $\Delta$ are thoroughly discussed.

5.1. *Practical implementation.* The computation of the three estimators $\hat{m}, \hat{v}, \hat{w}$ is obvious. However, the choice and computation of the two inverse functions $\Phi^{(+)}$ and $\Phi^{(-)}$, defined in Appendix G, need further explanations.

First, one has to compute $d^{(a)}$ defined in Equation (109). To avoid numerical errors, we replaced the condition $\kappa(m, w) = 4r_+r_-$ by the condition $|\kappa - 4r_+r_-| < 10^{-4}$. Moreover, we replaced $(4r_+r_-)^2 - 4r_+r_- + \kappa(m, w)$ by $\max\{0, (4r_+r_-)^2 - 4r_+r_- + \kappa(m, w)\}$ to ensure that the square-root stays in $\mathbb{R}_+$. Remark that in those two cases, we have $d^{(+)} = d^{(-)}$. In particular, there is no issue regarding the choice of $\Phi^{(+)}$ or $\Phi^{(-)}$. Out of those two cases, we use $\Phi^{(-)}$ whenever

$$ r_+ \ge 1/2 \quad \text{or} \quad \kappa(m, w) \ge 4r_+r_- \quad \text{or} \quad d^{(+)}(m, w) > 2r_- $$

The first two conditions are chosen according to Proposition G.3. The last one is supported by the fact that, if $r_+ < 1/2$ and $\kappa(m, w) < 4r_+r_-$ then $D(\lambda, p) < 2r_-$ (see Proposition G.1) and $d^{(-)}(m, w) < d^{(+)}(m, w)$. Finally, if all the conditions above are not satisfied then the choice between $a = +$ and $a = -$ is arbitrarily made.

Second, one has to compute $\phi_1^{(a)}$. On the one hand, to avoid numerical errors and the choice between $a = +$ and $a = -$, we replaced the condition $r_+ = 1/2$ by $|r_+ - r_-| < 10^{-3}$. On the other hand, when $\phi_1^{(a)}$ is applied to the estimators $\hat{m}$ and $\hat{w}$, one may end up with a negative value. Yet, $\phi_1^{(a)}$ should be non negative. Hence, we replaced $\phi_1^{(a)}(m, w)$ by $|\phi_1^{(a)}(m, w)|$ (for instance in Equation (113)).

Finally, when $\Phi^{(+)}$ or $\Phi^{(-)}$ are applied to the estimators $(\hat{m}, \hat{v}, \hat{w})$, one may end up with estimators $(\hat{\mu}, \hat{\lambda}, \hat{p})$ which do not belong to the admissible set $\Lambda$ defined above Equation (9). In that case we chose to clip the values. For instance, clip the value of $\lambda$ to 0 if negative, or to 1 if strictly larger than 1.

The numeric experiments were made using Julia programming language and a package should be available for the next revision of the paper.

**5.2. Choice of the tuning parameter.** From the bound obtained in Equation (12), the optimal choice of $\Delta$ is $\Delta = -\log(T)/(2\log(1-\lambda))$. In that case,

$$ \frac{(1-\lambda)^\Delta}{\Delta} + \sqrt{\frac{\Delta}{T}} = \frac{-2\log(1-\lambda)}{T^{1/2}\log(T)} + \frac{\sqrt{\log(T)}}{T^{1/2}}. $$

Since the value of $\lambda$ is not known, the choice $\Delta = \log(T)$ is made by default. However, the choice $\Delta = 1$ seems to give the best results among several tests and across several choices of parameters (see Figure 1). This is the reason why most of the plots of the next section are made with $\Delta = 1$.
<figure><img src="image_4.png" /><figcaption>FIG 1. Absolute estimation error for the six estimators and their theoretical limits. The y-axis is in log-scale. Each line or mark correspond to a median computed over $N_{simu} = 100$ simulations. The panels correspond to the choices $\Delta = \log(T)$ and $\Delta = 1$ from left to right.</figcaption></figure>

5.3. Numerical results. Let us first give the framework of our simulations. If not specified otherwise, the following values are used:

$$N = 500, \quad r_+ = .5, \quad \beta = .5, \quad \lambda = .5, \quad p = .5, \quad \Delta = 1, \quad N_{simu} = 1000.$$

Furthermore, the fractions of excitatory and inhibitory components are chosen as $r_+^N = \lceil r_+ N \rceil$ and $r_-^N = N - \lceil P_+ \rceil$.

In all the plots, the performance of one (or several) of our estimators is displayed. Let us denote $\hat{\vartheta}$ one of these estimators and $\vartheta$ its corresponding true parameter value. The solid lines correspond to the median of $|\hat{\vartheta} - \vartheta|$ computed over $N_{simu}$ simulations and plotted as a function of the time horizon $T$. Furthermore, we add marks (circles or horizontal bars) on the right end of the plot. These marks correspond to the median of $|\hat{\vartheta}_\infty - \vartheta|$ computed over $N_{simu}$ simulations, where $\hat{\vartheta}_\infty$ is the theoretical limit of $\hat{\vartheta}$ as $T \to \infty$ while the environment $\theta$ is fixed (remind Equation (31)). Of course, these quantities are unknown in practice but we are able to compute them here since the parameters used for the simulation are known. Finally, the corresponding limit estimators of $\mu$, $\lambda$ and $p$ are defined by:

$$\langle \hat{\mu}_{\infty}, \hat{\lambda}_{\infty}, \hat{p}_{\infty} \rangle = \Phi^{(a)} (\hat{m}_{\infty}, \hat{v}_{\infty}, \hat{w}_{\infty}),$$

where the choice of $a \in \{+, -\}$ is made according to Section 5.1.

Figure 1 compares the performance of all the estimators. The two panels correspond to two choices of $\Delta$ that will be discussed below. On the two panels, it is clear that the convergence of $\hat{m}$ is faster than all the others (which is expected from our analysis). Furthermore, for large $T$, $\hat{m}$, $\hat{v}$ and $\hat{\lambda}$ are really close to their theoretical limits $\hat{m}_\infty$, $\hat{v}_\infty$ and $\hat{\lambda}_\infty$ respectively. Note that this is not the case for $\hat{w}$, $\hat{\mu}$ and $\hat{p}$. In particular, it seems that $\hat{p}$ has the slowest convergence rate.

Figure 1 gives also a comparison of the performance of our estimators with respect to the tuning parameter $\Delta$. Either $\Delta$ is chosen as a function of $T$ by $\Delta = \log(T)$ or $\Delta$ is fixed to the value 1. Obviously, the estimators $\hat{m}$, $\hat{v}$ for finite $T$ and all the limit estimators $\hat{\vartheta}_\infty$ do not depend on $\Delta$ so the differences between the two plots are due to randomness only. However, it seems that the convergence of $\hat{w}$ is significantly faster when $\Delta = 1$. In turn, it implies the same improvement for the triplet $(\hat{\mu}, \hat{\lambda}, \hat{p})$.

Figure 2 gives an overview of the performance of the estimator $\hat{p}$ as one of the parameters $(N, r_+, \beta, \lambda$ or $p)$ varies. The choice to vary $\beta$ instead of $\mu = \lambda\beta$ is made because the set of admissible values of $\beta$ is independent of $\lambda$ (which is not the case for $\mu$). Remark that when $p = 0.5$, which is the case for most of the panels in Figure 2, the value 0.25
<figure><img src="image_5.png" /><figcaption>FIG 2. Absolute estimation error of $\hat{p}$ and its theoretical limit. Each line or mark corresponds to a median computed over $N_{simu} = 1000$ simulations. The panels correspond to different choices of varying parameter (the non-varying parameters are chosen according to the default values given in Section 5.3). The values of the varying parameter are given by the color legends.</figcaption></figure>

correspond to the median of the absolute error of the most naive estimator: estimate $p$ by a random uniform value in $[0, 1]$. Overall, the estimation error of the limit $\hat{p}_{\infty}$ is usually of the order of 0.01 and that of the estimator $\hat{p}$ computed for $T = 1000$ is approximately 20 times larger. As $N$ increases, the performance of $\hat{p}$ for $T$ fixed deteriorates. This phenomenon is encoded in the factor $\sqrt{N}/T$ appearing in Theorem 2.1 for example. Also, the marks corresponding to the estimation errors of the theoretical limit $\hat{p}_{\infty}$ for different values of $N$ are ordered as expected: the estimation error goes to 0 as $N$ goes to infinity. As $r_+$ varies, the performance of $\hat{p}$ seems not to vary too much. However, it is important to note that the case $r_+ = 0.4$ corresponds to a case where the choice of $a \in \{+, -\}$ is not obvious and so it is arbitrary (see Section 5.1). Half of the time, the wrong $a$ is chosen. In particular, it is the reason why the mark corresponding to the limit estimator is around 0.06 instead of 0.01 for all the other cases. As $\beta$ varies, the performances of $\hat{p}$ and $\hat{p}_{\infty}$ seem not to vary too much. This is quite expected. As $\lambda$ increases towards 1, the performances of $\hat{p}$ and $\hat{p}_{\infty}$ decreases. This is expected since $\lambda \to 1$ implies that the strength of interactions in the system
goes to 0. In particular, it is more and more difficult to estimate $p$ which is closely related to the interactions. Nevertheless, contrarily to our upper-bounds which goes to infinity as $\lambda \to 0$ (see the remark below Theorem 2.1), the performance is good for $\lambda = 0.1$. Finally, the performance of $\hat{p}$ as $p$ varies is more complex to analyze. For instance, for $p = 0.9$, our method gives $\hat{p} = 1$ most of the time whence the eventually constant violet curve at $y = 0.1$. In the case $p = 0.1$, the performance of $\hat{p}$ seems to be poor despite the fact that its limit $\hat{p}_{\infty}$ (blue mark) achieves good results. This phenomenon is also expected because as $p \to 0$ the interactions in the system vanishes.

**Acknowledgments.** This research has been conducted while J.C. was in Statify team at Centre Inria de l'Université Grenoble Alpes.

**Funding.** G.O. was supported by the Serrapilheira Institute (grant number Serra - 2211-42049), FAPERJ (grants E-26/201.397/2021 and E-26/204.532/2024) and CNPq (grants 303166/2022-3). J.C. and E.L. were supported by ANR-19-CE40-0024 (CHAlenges in MAthematical NEuroscience).

## REFERENCES

BASU, S. and MICHAILIDIS, G. (2015). Regularized estimation in sparse high-dimensional time series models. *The Annals of Statistics* **43** 1535–1567. https://doi.org/10.1214/15-AOS1315

BRESLER, G. (2015). Efficiently Learning Ising Models on Arbitrary Graphs. In *Proceedings of the Forty-Seventh Annual ACM Symposium on Theory of Computing. STOC '15* 771–782. Association for Computing Machinery, New York, NY, USA. https://doi.org/10.1145/2746539.2746631

BYOL KIM, S. L. and KOLAR, M. (2021). Two-Sample Inference for High-Dimensional Markov Networks. *Journal of the Royal Statistical Society Series B: Statistical Methodology* **83** 939–962.

CASELLA, G. and BERGER, R. (2024). *Statistical inference*. CRC press.

CHEVALLIER, J. and OST, G. (2024). Community detection for binary graphical models in high dimension. *arXiv preprint arXiv:2411.15627*.

COMETS, F., FERNÁNDEZ, R. and FERRARI, P. A. (2002). Processes with long memory: Regenerative construction and perfect simulation. *The Annals of Applied Probability* **12** 921–943. https://doi.org/10.1214/aoap/1031863175

DE SANTIS, E. and PICCIONI, M. (2015). One-Dimensional Infinite Memory Imitation Models with Noise. *Journal of Statistical Physics* **161** 346–364. https://doi.org/10.1007/s10955-015-1335-5

DELATTRE, S. and FOURNIER, N. (2016). Statistical inference versus mean field limit for Hawkes processes. *Electronic Journal of Statistics* **10** 1223–1295.

DUARTE, A., GALVES, A., LÖCHERBACH, E. and OST, G. (2019). Estimating the interaction graph of stochastic neural dynamics. *Bernoulli* **25** 771–792.

EICHLER, M. (2012). Graphical modelling of multivariate time series. *Probability Theory and Related Fields* **153** 233–268.

FERNÁNDEZ, R., FERRARI, P. A. and GALVES, A. (2001). Coupling, renewal and perfect simulation of chains of infinite order. Notes for a minicourse given at the Vth Brazilian School of Probability.

FERRARI, P. A. (1990). Ergodicity for Spin Systems with Stirrings. *The Annals of Probability* **18** 1523–1538. https://doi.org/10.1214/aop/1176990629

FERRARI, P. A., MAASS, A., MARTÍNEZ, S. and NEY, P. (2000). Cesàro mean distribution of group automata starting from measures with summable decay. *Ergodic Theory and Dynamical Systems* **20** 1657–1670. https://doi.org/10.1017/S0143385700000924

GAITONDE, J., MOITRA, A. and MOSSEL, E. (2024). Bypassing the Noisy Parity Barrier: Learning Higher-Order Markov Random Fields from Dynamics. *arXiv preprint arXiv:2409.05284*.

KOMLÓS, J., MAJOR, P. and TUSNÁDY, G. (1975). An Approximation of Partial Sums of Independent RV's and the Sample DF. I. *Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete* **32** 111–131. https://doi.org/10.1007/S0143385700000924

LAURITZEN, S. L. (1996). *Graphical Models*. Oxford University Press.

LEHMANN, E. L. and CASELLA, G. (1998). *Theory of Point Estimation*, 2nd ed. Springer Texts in Statistics. Springer. https://doi.org/10.1007/b98854

LERASLE, M. and TAKAHASHI, D. Y. (2016). Sharp oracle inequalities and slope heuristic for specification probabilities estimation in discrete random fields. *Bernoulli* **22** 325-344.
LIU, C. (2020). Statistical Inference for a partially observed interacting System of Hawkes Processes. *Stochastic Processes and Their Applications* **130** 5636-5694.

MONTANARI, A. and PEREIRA, J. (2009). Which graphical models are difficult to learn? In *Advances in Neural Information Processing Systems* (Y. BENGIO, D. SCHUURMANS, J. LAFFERTY, C. WILLIAMS and A. CU-LOTTA, eds.) **22**. Curran Associates, Inc.

RAVIKUMAR, P., WAINWRIGHT, M. J. and LAFFERTY, J. D. (2010). High-dimensional Ising model selection using $\ell_1$-regularized logistic regression. *Ann. Stat.* **38** 1287-1319. https://doi.org/10.1214/09-AOS691

REYNAUD-BOURET, P., RIVOIRARD, V. and TULEAU-MALOT, C. (2013). Inference of functional connectivity in Neurosciences via Hawkes processes. In *2013 IEEE Global Conference on Signal and Information Processing* 317-320.

STROGATZ, S. H. (2001). Exploring complex networks. *Nature* **410** 268-276.
APPENDIX A: PROOF OF THEOREM 3.3

The stationarity of the process X follows from the fact that its construction is invariant under time shift. Hence, to conclude the proof, it remains to show that, conditionally on each realization of the random environment $\theta$, X is Markovian and has the transition probabilities given by (2) and (1). The Markovianity of the process and the conditional independence of its coordinates follow from its representation as a function of an i.i.d. sequence. The compatibility of the transition probabilities follows from simple computations.

*Markovianity of the process.* Observe that (29) implies that for each $i \in \{1, \dots, N\}$, there exists a deterministic function $f_{\theta,i} : \{0,1\}^N \times \{0,1\} \times \{0,1, \dots, N\} \to \{0,1\}$ such that for all $t \in \mathbb{Z}$,

$$X_{i,t} = f_{\theta,i}(X_{t-1}, \xi_{i,t}, J_{i,t}).$$

Then defining the function

$$f_{\theta}(x, \xi, J) = (f_{\theta,1}(x, \xi_1, J_1), \dots, f_{\theta,N}(x, \xi_N, J_N)),$$

where $x \in \{0,1\}^N$, $\xi = (\xi_1, \dots, \xi_N) \in \{0,1\}^N$ and $J = (J_1, \dots, J_N) \in \{0,1, \dots, N\}^N$, we can write for all $t \in \mathbb{Z}$,

$$X_t = f_{\theta}(X_{t-1}, \xi_t, J_t),$$

where $\xi_t = (\xi_{1,t}, \dots, \xi_{N,t})$ and $J_t = (J_{1,t}, \dots, J_{N,t})$. This representation together with the fact the sequence $(\xi_t, J_t)_{t \in \mathbb{Z}}$ is i.i.d ensures that $(X_t)_{t \in \mathbb{Z}}$ is a Markov chain. Moreover, since the pairs $(\xi_{1,t}, J_{1,t})$, $\dots$, $(\xi_{N,t}, J_{N,t})$ are independent, the random variables $X_{i,t} = f_{\theta,i}(X_{t-1}, \xi_{i,t}, J_{i,t})$, $i \in \{1, \dots, N\}$, are clearly conditionally independent given $X_{t-1}$.

*Transition probability.* It remains to show that $\mathbb{P}_{\theta}(X_{i,t} = 1 | X_{t-1} = x) = p_{\theta,i}(x)$ for all $x \in \{0,1\}^N$, where $p_{\theta,i}(x)$ is given by (1). To see that, write $z = (i,t)$ and observe that (29) allows us to write,

$$\mathbb{P}_{\theta}(X_z = 1 | X_{t-1} = x) = \mathbb{P}_{\theta}(J_z = 0, \xi_z = 1) + \sum_{j \in \mathcal{P}_{+}} \theta_{ij} x_j \mathbb{P}_{\theta}(J_z = j | X_{t-1} = x) \\ + \sum_{j \in \mathcal{P}_{-}} \theta_{ij} (1 - x_j) \mathbb{P}_{\theta}(J_z = j | X_{t-1} = x).$$

Using that $X_{t-1}$ is measurable with respect to $\sigma(\xi_{k,s}, J_{k,s}, s \le t-1, k \in \{1, \dots, N\})$ and the fact that both $J_z$ and $\xi_z$ are independent of all random variables $J_{k,s}$, $\xi_{k,s}$ for $s \le t-1$ and $k \in \{1, \dots, N\}$, we have that

$$\mathbb{P}_{\theta}(J_z = 0, \xi_z = 1 | X_{t-1} = x) = \mathbb{P}_{\theta}(J_z = 0, \xi_z = 1),$$

and also that for each $j \in \{1, \dots, N\}$,

$$\mathbb{P}_{\theta}(J_z = j | X_{t-1} = x) = \mathbb{P}_{\theta}(J_z = j) = \frac{(1 - \lambda)}{N}.$$

As a consequence, it follows that

$$\mathbb{P}_{\theta}(X = 1 | X_{t-1} = x) = \mathbb{P}_{\theta}(J_z = 0, \xi_z = 1) + \frac{(1 - \lambda)}{N} \sum_{j \in \mathcal{P}_{+}} \theta_{ij} x_j + \sum_{j \in \mathcal{P}_{-}} \theta_{ij} (1 - x_j).$$

Since the random variables $J_z$ and $\xi_z$ are also independent, and $\mathbb{P}_{\theta}(J_z = 0) = \lambda$, $\mathbb{P}_{\theta}(\xi_z = 1) = \beta$ and $\mu = \lambda\beta$, it follows that $\mathbb{P}_{\theta}(X_{i,t} = 1 | X_{t-1} = x) = p_{\theta,i}(x)$, concluding the proof.
APPENDIX B: COALESCENCE COUPLINGS

This section is devoted to the proof of two lemmas which relate independence and coalescence. It is structured as follows. First, the two lemmas are stated. Next, the construction of the coupling is detailed. After that, some properties of the coupling are provided. Finally, the proofs of the two lemmas are given.

**B.1. Two main lemmas.** Here are two lemmas that relate coalescence with independence via couplings. The tilde and hat versions of $X$ defined below respect some additional spatial independence properties. The objective of the two lemmas is to prove that these versions can be coupled in such a way that they are equal except on an event related with coalescence (which in turn we expect to be of small probability).

LEMMA B.1 (All independent coupling). Let $k \ge 2$ and consider $k$ different points $z_1, \dots, z_k$ in $\mathcal{Z}$. There exist random variables $\tilde{X}_{z_1}, \dots, \tilde{X}_{z_k}$ satisfying the following properties: for each $i \in \{1, \dots, k\}$,

(i) $\tilde{X}_{z_i} \stackrel{d}{=} X_{z_i}$ (i.e., $\tilde{X}_{z_i}$ and $X_{z_i}$ have the same distribution);

(ii) $\tilde{X}_{z_i}$ is independent of $(X_{z_j}, \tilde{X}_{z_j})_{j \in \{1, \dots, k\} \setminus \{i\}}$;

(iii) $\{\tilde{X}_{z_i} \neq X_{z_i}\} \subset \cup_{j \in \{1, \dots, k\} \setminus \{i\}} \{z_i \rightsquigarrow z_j\}$.

LEMMA B.2 (Blockwise independence). Let $z_1, \dots, z_4$ be four different points in $\mathcal{Z}$. There exist random variables $\tilde{X}_{z_1}, \dots, \tilde{X}_{z_4}$ and $\hat{X}_{z_1}, \hat{X}_{z_2}$ satisfying the following properties:

(i) $(\hat{X}_{z_1}, \hat{X}_{z_2}) \stackrel{d}{=} (X_{z_1}, X_{z_2})$;

(ii) $(\hat{X}_{z_1}, \hat{X}_{z_2})$ is independent of $(X_{z_3}, X_{z_4}, \tilde{X}_{z_3}, \tilde{X}_{z_4})$;

(iii) $\{\hat{X}_{z_1} \neq X_{z_1}\} \subset \cup_{k=3}^4 \{z_1 \rightsquigarrow z_k\} \cup (\{\tau_{z_1}^R < t_3 \lor t_4\} \cap (\cup_{k=3}^4 \{z_2 \rightsquigarrow z_k\}))$,

(iv) $\{\hat{X}_{z_2} \neq X_{z_2}\} \subset \cup_{k=3}^4 \{z_2 \rightsquigarrow z_k\} \cup (\{\tau_{z_2}^R < t_3 \lor t_4\} \cap (\cup_{k=3}^4 \{z_1 \rightsquigarrow z_k\}))$.

(v) $\hat{X}_{z_1}$ is independent of $\tilde{X}_{z_2}$ and $\hat{X}_{z_2}$ is independent of $\tilde{X}_{z_1}$.

Moreover, for each $i \in \{1, \dots, 4\}$,

(iv) $\tilde{X}_{z_i} \stackrel{d}{=} X_{z_i}$;

(v) $\tilde{X}_{z_i}$ is independent of $(\tilde{X}_{z_j}, X_{z_j}, I^{z_j})_{j \in \{1, \dots, 4\} \setminus \{i\}}$;

(vi) $\{\tilde{X}_{z_i} \neq X_{z_i}\} \subset \cup_{j \in \{1, \dots, 4\} \setminus \{i\}} \{z_i \rightsquigarrow z_j\}$.

REMARK B.3. Let us describe with some words Item (iii) above: if the hat version $\hat{X}_{z_1}$ is different from the standard version $X_{z_1}$, then we know that

• $z_1$ coalesces with $z_3$ or $z_4$, or

• $z_2$ coalesces with $z_3$ or $z_4$ and the (backward) regeneration time of $z_1$ is less than the starting times of $z_3$ or $z_4$.

**B.2. Construction of the couplings.** Let us denote $(J, \xi)$ the sequence $(J_z, \xi_z)_z$ defined in the backward regeneration representation of the process $X$. Since we are working for a fixed matrix $\theta$ in this section, we omit its dependence in what follows. Due to the backward regeneration representation there exist measurable functions $f, g, r$ and $\tau$ (which depend on $\theta$) such that, for all $z \in \mathcal{Z}$,

$$I^z = (I_s^z)_{s \in \mathbb{Z}} = f(z, J) \text{ and } X_z = g(I^z, \xi_{r(I^z)}),$$
where $r(I^z) = (I_{\tau(I^z)}^z, \tau(I^z))$ is the regenerating site associated with $z$ and $\tau(I^z) = \inf\{s \in \mathbb{Z} : I_s^z > 0\}$ is the last time that the random walk $I^z$ is not in the state 0. In what follows, we denote $I = (I^z)_{z \in \mathbb{Z}}$.

Now, consider 5 independent sequences

$$(J^{(1)}, \xi^{(1)}), \dots, (J^{(5)}, \xi^{(5)}),$$

all distributed as the sequence $(J, \xi)$, and which are furthermore independent of $(J, \xi)$. By convenience of notation, let us denote $(J^{(0)}, \xi^{(0)}) := (J, \xi)$. In other words, the layer (0) corresponds to the original process, and we add five i.i.d. layers on top of that. We then define the backward random walks $I^{(i)} = (I^{(i),z})_{z \in \mathbb{Z}}$ by $I^{(i),z} = f(z, J^{(i)})$, for each $0 \le i \le 5$. Clearly, the collections of random walks $I^{(0)}, \dots, I^{(5)}$ are independent. We will need the layers (0) – (4) to construct the tilde versions, and the layers (0) and (5) to construct the hat versions.

*Construction of the tilde versions.* For $i \in \{1, \dots, 4\}$, let us denote

$$\tau_i^c = \sup \left\{ \tau_{z_i, z_j}^c : j \in \{1, \dots, 4\} \setminus \{i\} \right\},$$

the first time at which $I^{z_i}$ coalesces with (at least) one of the 3 other random walks (observe that $\tau_i^c = -\infty$ if $I^{z_i}$ does not coalesce with any of the other 3 random walks), and denote $z_i^c = (I_{\tau_i^c}^{z_i}, \tau_i^c)$ the site of coalescence in case $\tau_i^c > -\infty$. Then, for all $1 \le i \le 4$ and $s \in \mathbb{Z}$, we define

$$(39) \qquad \tilde{I}_s^{z_i} = I_s^{(0),z_i} \mathbf{1}_{s > \tau_i^c} + I_s^{(i),z_i} \mathbf{1}_{s \le \tau_i^c}.$$

In words, the random walk $\tilde{I}^{z_i}$ follows the random walk $I^{z_i}$ until it coalesces with one of the random walks $I^{z_j}, j \in \{1, \dots, 4\} \setminus \{i\}$. After this time, the random walk $\tilde{I}^{z_i}$ follows the independent random walk $I^{(i),z_i^c}$ associated with the site of coalescence $z_i^c$. With the notation $\tilde{I}_s^{z_i} = i \cdot 1_{s \le \tau_i^c} \in \{0, \dots, 4\}$, Equation (39) rewrites as

$$\tilde{I}_s^{z_i} = I_s^{(0),z_i} \mathbf{1}_{\tilde{I}_s^{z_i}=0} + I_s^{(i),z_i} \mathbf{1}_{\tilde{I}_s^{z_i}=i},$$

and some computations give the following easy recursion formula:

$$(40) \qquad \tilde{I}_{s-1}^{z_i} = J_{(\tilde{I}_s^{z_i}, s)}^{(\tilde{I}_s^{z_i})}, \text{ for all } s \le t_i,$$

which can be compared with Equation (28). Moreover, remark that the layer processes $\tilde{L}$ are measurable functions of the state processes $\tilde{I}$ since

$$\tau_{z_i, z_j}^c = \sup \{s \in \mathbb{Z} : I_s^{z_i} = I_s^{z_j} \notin \{0, \infty\}\} = \sup \{s \in \mathbb{Z} : \tilde{I}_s^{z_i} = \tilde{I}_s^{z_j} \notin \{0, \infty\}\}.$$

*Construction of the hat versions.* Let us denote, for $i = 1, 2$,

$$\hat{\tau}_i = \sup \left\{ \tau_{z_i, z_j}^c : j \in \{3, 4\} \right\},$$

and denote $\hat{z}_i^c = (I_{\hat{\tau}_i}^{z_i}, \hat{\tau}_i)$ the corresponding site of coalescence in case $\hat{\tau}_i > -\infty$. For instance, $\hat{\tau}_1$ corresponds to the first time at which $I^{z_1}$ coalesces with either $I^{z_3}$ or $I^{z_4}$. Then, for $i = 1, 2$, we define the auxiliary process, for all $s \in \mathbb{Z}$,

$$\hat{I}_s^{(5),z_i} = I_s^{z_i} \mathbf{1}_{s > \hat{\tau}_i} + I_s^{(5),\hat{z}_i} \mathbf{1}_{s \le \hat{\tau}_i}.$$

The construction is not yet done since these two random walks may be missing some coalescence in order to mimic the joint distribution of $(X_{z_1}, X_{z_2})$. To circumvent this problem, we define

$$\hat{\tau}_{1,2} = \sup \left\{ s \in \mathbb{Z} : \hat{I}_s^{(5),z_1} = \hat{I}_s^{(5),z_2} \notin \{0, \infty\} \right\},$$
and

$$ \hat{I}^{z_i} = \hat{I}^{(5),z_i}, \text{ if } \hat{\tau}_{1,2} \le \hat{\tau}_i, $$

while we put, for all $s \in \mathbb{Z}$,

$$ (41) \qquad \hat{I}_s^{z_i} = I_s^{z_i} \mathbf{1}_{s>\hat{\tau}_{1,2}} + \hat{I}_s^{(5),z_{3-i}} \mathbf{1}_{s\le\hat{\tau}_{1,2}}, \text{ if } \hat{\tau}_{1,2} > \hat{\tau}_i. $$

For instance, if $I^{z_1}$ first coalesces, say, with $I^{z_3}$, then it switches to its hat-version $\hat{I}^{z_1}$ and remains stuck to this version forever. However, if $I^{z_2}$ hits $\hat{I}^{z_1}$ before hitting $I^{z_3}$ or $I^{z_4}$, then it coalesces with $\hat{I}^{z_1}$ and remains stuck to it forever. In words, we do not modify the random walk with the first (the largest in time) coalescence time $\hat{\tau}_i$, but we possibly modify the other one in between the times $\hat{\tau}_1$ and $\hat{\tau}_2$. Like Equation (40), we have the following recursion:

$$ (42) \qquad \hat{I}_{s-1}^{z_i} = J_{(\hat{I}_s^{z_i}, s)}^{(\hat{I}_s^{z_i})}, \text{ for all } s \le t_i, $$

where

$$ (43) \qquad \hat{L}_s^{z_1} = 5 \cdot 1_{s \le \hat{\tau}_1} + 5 \cdot 1_{\hat{\tau}_1 < s \le \hat{\tau}_{1,2}} \mathbf{1}_{\hat{\tau}_{1,2} \le \hat{\tau}_2} $$

with a similar definition for $\hat{L}_s^{z_2}$. In other words, in any case, process 1 is in layer 5 starting from time $\hat{\tau}_1$, that is $\hat{L}_s^{z_1} = 5$ for $s \le \hat{\tau}_1$. But, if process 1 coalesces with process 2 (at time $\hat{\tau}_{1,2}$) and if process 2 is already evolving on layer 5 at that time, then process 1 switches to layer 5 at that coalescence time $\hat{\tau}_{1,2}$.

Once again, the layer processes $\hat{L}^{z_1}$ and $\hat{L}^{z_2}$ are measurable functions of the state processes $(\hat{I}^{z_1}, \hat{I}^{z_2}, I^{z_3}, I^{z_4})$.

**B.3. Properties of the couplings.** First, let us formalize the fact that the processes $I$ do not depend on the whole sequences $J^{(0)}, \dots, J^{(5)}$ but merely on a small subset of those. Let $\gamma \in \{0, 1, \dots, N, \infty\}^{\mathbb{Z}}$ and $\ell \in \{0, \dots, 5\}^{\mathbb{Z}}$ be some generic trajectories of the processes $I$ and $L$. Let us then define $J_{\gamma}^{\ell} = (J_{(\gamma_s, s)}^{(\ell_s)})_{s \in \mathbb{Z}} \in \{0, 1, \dots, N\}^{\mathbb{Z}}$ and, for all $z = (i, t) \in \mathbb{Z}$, $f_z^{\ell} : \{0, 1, \dots, N\}^{\mathbb{Z}} \to \{0, 1, \dots, N, \infty\}^{\mathbb{Z}}$ given by

$$ (f_z^{\ell}(J_{\gamma}^{\ell}))_s = \begin{cases} \infty & \text{if } s > t \\ i & \text{if } s = t \\ J_{(\gamma_{s+1}, s+1)}^{(\ell_{s+1})} & \text{else.} \end{cases} $$

Finally, let us remark that:

1. the definition of $I^{z_i}$ implies that $\{I^{z_i} = \gamma\} = \{f_{z_i}^{\ell}(J_{\gamma}^{\ell}) = \gamma\}$ where $\ell_s^0 = 0$ for all $s \in \mathbb{Z}$;

2. Equation (40) and the measurability of the layer processes imply that for all $\tilde{\gamma}^1, \dots, \tilde{\gamma}^4 \in \{0, 1, \dots, N, \infty\}^{\mathbb{Z}}$, there exist $\tilde{\ell}^1, \dots, \tilde{\ell}^4 \in \{0, \dots, 4\}^{\mathbb{Z}}$ such that

$$ \cap_{i=1}^{4} \{\tilde{I}^{z_i} = \tilde{\gamma}^i\} = \cap_{i=1}^{4} \{\tilde{I}^{z_i} = \tilde{\gamma}^i, \tilde{L}^{z_i} = \tilde{\ell}^i\} = \cap_{i=1}^{4} \{f_{z_i}^{\ell}(J_{\tilde{\gamma}^i}^{\ell}) = \tilde{\gamma}^i\}; $$

3. Equation (42) and the measurability of the layer processes imply that for all choices $\hat{\gamma}^1, \hat{\gamma}^2, \gamma^3, \gamma^4, \tilde{\gamma}^3, \tilde{\gamma}^4 \in \{0, 1, \dots, N, \infty\}^{\mathbb{Z}}$, there exist $\hat{\ell}^1, \hat{\ell}^2 \in \{0, 5\}^{\mathbb{Z}}$ and $\tilde{\ell}^3, \tilde{\ell}^4 \in \{0, 3, 4\}^{\mathbb{Z}}$ such that

$$ \begin{align*} A &:= (\cap_{i=1}^{2} \{\hat{I}^{z_i} = \hat{\gamma}^i\}) \cap (\cap_{j=3}^{4} \{I^{z_j} = \gamma^j, \tilde{I}^{z_j} = \tilde{\gamma}^j\}) \\ &= (\cap_{i=1}^{2} \{\hat{I}^{z_i} = \hat{\gamma}^i, \hat{L}^{z_i} = \hat{\ell}^i\}) \cap (\cap_{j=3}^{4} \{I^{z_j} = \gamma^j, \tilde{I}^{z_j} = \tilde{\gamma}^j, \tilde{L}^{z_j} = \tilde{\ell}^j\}) \\ &= (\cap_{i=1}^{2} \{f_{z_i}^{\ell}(J_{\hat{\gamma}^i}^{\ell}) = \hat{\gamma}^i\}) \cap (\cap_{j=3}^{4} \{f_{z_j}^{\ell}(J_{\hat{\gamma}^j}^{\ell}) = \gamma^j, f_{z_j}^{\ell}(J_{\tilde{\gamma}^j}^{\ell}) = \tilde{\gamma}^j\}). \end{align*} $$
First, let us summarize the coupling properties of the processes I.

**PROPOSITION B.4.** The processes I, $\tilde{I}$ and $\hat{I}$ satisfy the following coupling properties:

(i) for all $i = 1, \dots, 4$, $\{\tilde{I}^{z_i} \neq I^{z_i}\} = \cup_{j \in \{1, \dots, 4\} \setminus \{i\}} \{z_i \rightsquigarrow z_j\}$,

(ii) we have that

$$ \{\hat{I}^{z_1} \neq I^{z_1}\} \subset \bigcup_{k=3}^{4} \{z_1 \rightsquigarrow z_k\} \cup \left( \bigcup_{k=3}^{4} \{z_2 \rightsquigarrow z_k\} \cap \{\tau_{z_1}^{R} < t_3 \lor t_4\} \right) $$

and

$$ \{\hat{I}^{z_2} \neq I^{z_2}\} \subset \bigcup_{k=3}^{4} \{z_2 \rightsquigarrow z_k\} \cup \left( \bigcup_{k=3}^{4} \{z_1 \rightsquigarrow z_k\} \cap \{\tau_{z_2}^{R} < t_3 \lor t_4\} \right). $$

**PROOF.** This follows directly from the construction of the processes and the remark that

$$ \{\hat{\tau}_i < \hat{\tau}_{1,2} \le \hat{\tau}_{3-i}\} \subset \bigcup_{k=3}^{4} \{z_{3-i} \rightsquigarrow z_k\} \cap \{\tau_{z_i}^{R} < t_3 \lor t_4\}. $$

<figure><img src="image_6.png" /></figure>

We are now in position to prove the following independence properties of the processes I.

**PROPOSITION B.5.** The processes I, $\tilde{I}$ and $\hat{I}$ satisfy the following independence properties:

(i) for all $i$, $\tilde{I}^{z_i}$ is independent of $(I^{z_j}, \tilde{I}^{z_j})_{j \neq i}$;

(ii) $(\hat{I}^{z_1}, \hat{I}^{z_2})$ is independent of $(I^{z_j}, \tilde{I}^{z_j})_{j=3,4}$.

(iii) $\hat{I}^{z_1}$ is independent of $\tilde{I}^{z_2}$, and $\hat{I}^{z_2}$ is independent of $\tilde{I}^{z_1}$.

**PROOF.**

*Proof of (i).* Assume without loss of generality that $i = 1$. Let $\tilde{\gamma}^1, \dots, \tilde{\gamma}^4, \gamma^2, \dots, \gamma^4 \in \{0, 1, \dots, 2N, \infty\}^{\mathbb{Z}}$, and let $\tilde{\ell}_1, \dots, \tilde{\ell}_4$ be the associated layer processes, which are uniquely defined being measurable functions of $\tilde{\gamma}^1, \dots, \tilde{\gamma}^4, \gamma^2, \dots, \gamma^4$. By the point 2 above, we have

$$ \begin{aligned} A &= (\cap_{i=1,...,4} \{\tilde{I}^{z_1} = \tilde{\gamma}^1\}) \cap (\cap_{j=2,...,4} \{I^{z_j} = \gamma^j\}) \\ &= \{f_{z_1}^{\tilde{\ell}_1}(J_{\tilde{\gamma}^1}^{\tilde{\ell}_1}) = \tilde{\gamma}^1\} \cap (\cap_{j=2}^{4} \{f_{z_j}^{\tilde{\ell}_j}(J_{\tilde{\gamma}^j}^{\tilde{\ell}_j}) = \tilde{\gamma}^j, f_{z_j}^{\ell_0}(J_{\tilde{\gamma}^j}^{\ell_0}) = \gamma^j\}) \end{aligned} $$

By construction, if $\mathcal{IP}(A)$ is non null, then the set $\{(\tilde{\gamma}_s^1, \tilde{\ell}_s^1, s), s \in \mathbb{Z}, \tilde{\gamma}_s^1 \notin \{0, \infty\}\}$ is disjoint from the sets $\{(\tilde{\gamma}_s^j, \tilde{\ell}_s^j, s), s \in \mathbb{Z}, \tilde{\gamma}_s^j \notin \{0, \infty\}\}$ and $\{(\gamma_s^j, 0, s), s \in \mathbb{Z}, \gamma_s^j \notin \{0, \infty\}\}$ for $j > 1$ so that the two events in the final expression of A are independent, and Item (i) follows.

*Proof of (ii).* Let $\hat{\gamma}^1, \hat{\gamma}^2, \gamma^3, \gamma^4, \tilde{\gamma}^3, \tilde{\gamma}^4 \in \{0, 1, \dots, N, \infty\}^{\mathbb{Z}}$, and let $\hat{l}^1, \hat{l}^2, \hat{l}^3, \hat{l}^4$ be the associated layer processes. By the point 3 above, we have

$$ \begin{aligned} A &= (\cap_{i=1}^{2} \{\hat{I}^{z_i} = \hat{\gamma}^i\}) \cap (\cap_{j=3}^{4} \{I^{z_j} = \gamma^j, \tilde{I}^{z_j} = \tilde{\gamma}^j\}) \\ &= (\cap_{i=1}^{2} \{f_{z_i}^{\tilde{\ell}_i}(J_{\tilde{\gamma}^i}^{\tilde{\ell}_i}) = \hat{\gamma}^i\}) \cap (\cap_{j=3}^{4} \{f_{z_j}^{\ell_0}(J_{\tilde{\gamma}^j}^{\ell_0}) = \gamma^j, f_{z_j}^{\tilde{\ell}_j}(J_{\tilde{\gamma}^j}^{\tilde{\ell}_j}) = \tilde{\gamma}^j\}). \end{aligned} $$
By construction, if $\mathcal{P}(A)$ is non null, then the sets $\{(\hat{\gamma}_s^1, \hat{l}_s^1, s), s \in \mathbb{Z}, \hat{\gamma}_s^1 \notin \{0, \infty\}\}$ and $\{(\hat{\gamma}_s^2, \hat{l}_s^2, s), s \in \mathbb{Z}, \hat{\gamma}_s^2 \notin \{0, \infty\}\}$ are disjoint from the sets $\{(\tilde{\gamma}_s^j, \tilde{l}_s^j, s), s \in \mathbb{Z}, \tilde{\gamma}_s^j \notin \{0, \infty\}\}$ and $\{(\gamma_s^j, 0, s), s \in \mathbb{Z}, \gamma_s^j \notin \{0, \infty\}\}$ for $j = 3, 4$ so that the two events in the final expression of A are independent, and Item (ii) follows.

*Proof of (iii).* Let $\hat{\gamma}^1, \tilde{\gamma}^2 \in \{0, 1, \dots, N, \infty\}^{\mathbb{Z}}$, and let $\hat{l}^1 \in \{0, 5\}^{\mathbb{Z}}$ and $\tilde{l}^2 \in \{0, 2\}^{\mathbb{Z}}$ be two fixed layer processes that are compatible with the event $\{\hat{I}^{z_1} = \hat{\gamma}^1\} \cap \{\tilde{I}^{z_2} = \tilde{\gamma}^2\}$. By compatible, we mean for example that supposing that $\hat{\gamma}^1$ and $\tilde{\gamma}^2$ meet at some time $s \le t_1 \lor t_2$, it is not possible to have $\hat{l}_u^1 = \tilde{l}_u^2$ for $u \le s$.

Then, analogously to the proof of point (ii) above, the sets $\{(\hat{\gamma}_s^1, \hat{l}_s^1, s), s \in \mathbb{Z}, \hat{\gamma}_s^1 \notin \{0, \infty\}\}$ and $\{(\tilde{\gamma}_s^2, \tilde{l}_s^2, s), s \in \mathbb{Z}, \tilde{\gamma}_s^2 \notin \{0, \infty\}\}$ are disjoint sets implying the independence. $\square$

**PROPOSITION B.6.** For all $i$, $\tilde{I}^{z_i}$ has the same distribution as $I^{z_i}$ and $(\hat{I}^{z_1}, \hat{I}^{z_2})$ has the same distribution as $(I^{z_1}, I^{z_2})$.

**PROOF.** Remind that the layers $J^{(0)}, \dots, J^{(5)}$ are i.i.d. so that it is easy to prove that $\tilde{I}^{z_i}$ is a backward random walk with the same transitions as $I^{z_i}$ (compare Equation (40) with Equation (28)) which in turn implies that they share the same distribution.

We now turn to the second part of the proof. For the marginals, the same argument as above applies (thanks to Equation (42)). Hence, $\hat{I}^{z_i}$ has the same distribution as $I^{z_i}$. Then, we show that the process $(\hat{I}_s^{z_1}, \hat{I}_s^{z_2})_s$ is a two-dimensional Markov chain (backwards in time) which has the same transitions as $(I_s^{z_1}, I_s^{z_2})_s$. Since both chains start from the same initial conditions, this implies the desired result.

Since transitions between $t_1 \lor t_2$ and $t_1 \land t_2$ only concern one of the two processes, evolving according to the right marginals, we do only need to consider transitions $s \to s-1$ for $s \le t_1 \land t_2$. Fix $i_1, i_2, j_1, j_2 \in \{0, \dots, N\}$. Let us first discuss the case $i_1 \neq i_2$.

Then for any $l_1, l_2 \in \{0, 5\}$, writing

$$A := \{\hat{I}_s^{z_1} = i_1, \hat{I}_s^{z_2} = i_2, \hat{L}_s^{z_1} = l_1, \hat{L}_s^{z_2} = l_2\},$$

we have that

$$\mathcal{P}(\hat{I}_{s-1}^{z_1} = j_1, \hat{I}_{s-1}^{z_2} = j_2 | A) = \mathcal{P}(J_{(i_1,s)}^{(l_1)} = j_1, J_{(i_2,s)}^{(l_2)} = j_2 | A).$$

By construction, $J_{(i_1,s)}^{(l_1)}$ and $J_{(i_2,s)}^{(l_2)}$ are independent (since $i_1 \neq i_2$), they have the same joint distribution as $(J_{(i_1,s)}^{(0)}, J_{(i_2,s)}^{(0)})$, and they are independent of A (since A depends only on decisions strictly after time s). Thus

$$\mathcal{P}(\hat{I}_{s-1}^{z_1} = j_1, \hat{I}_{s-1}^{z_2} = j_2 | A) = \mathcal{P}(J_{(i_1,s)}^{(0)} = j_1, J_{(i_2,s)}^{(0)} = j_2).$$

Summing over all possible choices of $l_1$ and $l_2$, this implies that

$$(44) \quad \mathcal{P}(\hat{I}_{s-1}^{z_1} = j_1, \hat{I}_{s-1}^{z_2} = j_2 | \hat{I}_s^{z_1} = i_1, \hat{I}_s^{z_2} = i_2) = \mathcal{P}(J_{(i_1,s)}^{(0)} = j_1, J_{(i_2,s)}^{(0)} = j_2).$$

We now discuss the case $i_1 = i_2$. Then necessarily $l_1 = l_2$ and thus, by our coalescence construction, $j_1 = j_2$. With the same notation for the set A as above and using the same independence argument, it still holds that

$$\mathcal{P}(\hat{I}_{s-1}^{z_1} = j_1 = \hat{I}_{s-1}^{z_2} | A) = \mathcal{P}(J_{(i_1,s)}^{(l_1)} = j_1 | A) = \mathcal{P}(J_{(i_1,s)}^{(l_1)} = j_1) = \mathcal{P}(J_{(i_1,s)}^{(0)} = j_1).$$

Summing over all possible values of $l_1$ implies that

$$(45) \quad \mathcal{P}(\hat{I}_{s-1}^{z_1} = j_1 = \hat{I}_{s-1}^{z_2} | \hat{I}_s^{z_1} = i_1 = \hat{I}_s^{z_2}) = \mathcal{P}(J_{(i_1,s)}^{(0)} = j_1).$$

The transition probabilities given by Equations (44) and (45) correspond exactly to the transition probabilities of $(I_s^{z_1}, I_s^{z_2})$ - see Equation (28). $\square$
## B.4. Proof of the two lemmas.

**PROOF OF LEMMA B.1.** Four tilde processes are constructed above so we write the proof for $k = 4$. Nevertheless, the proof can easily be generalized to any $k \ge 2$.

Let $z_1, \dots, z_4$ be four different sites in $\mathcal{Z}$. For all $i = 1, \dots, 4$, let $\tilde{I}^{z_i}$ be the backward random walks defined by (39), and $\tilde{L}^{z_i}$ be the associated layer process. We denote $\tilde{L}_{-\infty}^{z_i} = \lim_{s \to -\infty} \tilde{L}_s^{z_i} \in \{0, i\}$ the terminal layer. Then, remind the functions $g$ and $r$ defined in the beginning of Appendix B.2 and define, for all $i = 1, \dots, 4$,

$$ \tilde{X}_{z_i} = g\left(\tilde{I}^{z_i}, \xi_{r(\tilde{I}^{z_i})}^{(0)}\right) \mathbf{1}_{\tilde{L}_{-\infty}^{z_i}=0} + g\left(\tilde{I}^{z_i}, \xi_{r(\tilde{I}^{z_i})}^{(i)}\right) \mathbf{1}_{\tilde{L}_{-\infty}^{z_i}=i}. $$

In comparison, remind that $X_{z_i} = g(I^{z_i}, \xi_{r(I^{z_i})}^{(0)})$. Using the fact that the $\xi^{(j)}$'s are i.i.d. with Propositions B.5 and B.6, one deduces Items (i) and (ii). Finally, Item (iii) follows from the fact that $\{\tilde{X}^{z_i} \neq X^{z_i}\} \subset \{\tilde{I}^{z_i} \neq I^{z_i}\}$ and Proposition B.4. $\square$

**PROOF OF LEMMA B.2.** Let $z_1, \dots, z_4$ be four different sites in $\mathcal{Z}$. For all $i = 1, \dots, 4$, let us define $\tilde{X}_{z_i}$ as in the proof of Lemma B.1. In particular, Items (iv) to (vi) follow from Lemma B.1.

Then, for $i = 1, 2$, let $\hat{I}^{z_i}$ be the backward random walks defined by (41), and $\hat{L}^{z_i}$ be the associated layer process. We denote $\hat{L}_{-\infty}^{z_i} = \lim_{s \to -\infty} \hat{L}_s^{z_i} \in \{0, 5\}$ the terminal layer. Then, remind the functions $g$ and $r$ defined in the beginning of Appendix B.2 and define, for all $i = 1, \dots, 4$,

$$ \hat{X}_{z_i} = g\left(\hat{I}^{z_i}, \xi_{r(\hat{I}^{z_i})}^{(0)}\right) \mathbf{1}_{\hat{L}_{-\infty}^{z_i}=0} + g\left(\hat{I}^{z_i}, \xi_{r(\hat{I}^{z_i})}^{(5)}\right) \mathbf{1}_{\hat{L}_{-\infty}^{z_i}=5}. $$

Using the fact that the $\xi^{(j)}$'s are i.i.d., Item (i) follows from Proposition B.6, Items (ii) and (v) follow respectively from Items (ii) and (iii) of Proposition B.5. Finally, Items (ii) and (iii) follow from Item (ii) of Proposition B.4. $\square$

## APPENDIX C: COALESCENCE OF TWO OR MORE BACKWARD RANDOM WALKS

The aim of this section is to prove Proposition 3.5. Throughout the proof we will consider partitions of the sets of cardinal 2, 3 and 4. For a finite set $E$, we denote by $\mathcal{P}(E)$ the set of its partitions. For ease of notation, we consider a notation which we exemplify in the case of a set $\{a, b, c\}$ with three different elements. In this case, the set $\mathcal{P}(\{a, b, c\})$ has 5 elements which are:

* $\{\{a\}, \{b\}, \{c\}\}$ written as $a^\dagger b^\dagger c$,

* $\{\{a, b\}, \{c\}\}$, written as $ab^\dagger c$,

* $\{\{a, c\}, \{b\}\}$, written as $ac^\dagger b$,

* $\{\{a\}, \{b, c\}\}$, written as $a^\dagger bc$,

* $\{\{a, b, c\}\}$, written as $abc$.

The proof relies on the Markovian property of the backward random walks $(I^z)_{z \in \mathcal{Z}}$. From these backward random walks (BRW), we define a “status” random walk - such as $S_t^1$ below - which informs us whether the BRW has started, is currently running, or has already reached the cemetery. This, in turn, allows us to define a partition valued Markov chain, denoted $P_t$. To keep concise notation for the states of this partition valued Markov chain, we use the notation introduced above, which informs us on the coalescing status of the BRWs. The proof is then based on the computation of the transition probabilities of the partition Markov chain (some transition graphs are drawn for pedagogical purposes). Finally, the proof is concluded by the computation of the coalescence probabilities.

**PROOF OF PROPOSITION 3.5.** Without loss of generality, we assume in this proof that $t_1 \ge t_2 \ge t_3 \ge t_4$.
*Proof of Item (i).* Consider the backward random walks $I^{z_1}$ and $I^{z_2}$ from time $t = t_1$ to time $t = -\infty$. Let us introduce the following notation: in what follows, the letter 'S' stands for the state of a random walk; $S_t^1$ will be the state at time $t$ of the random walk associated to $z_1$, and $S_t^2$ the state at time $t$ of the random walk associated to $z_2$. We define $S_t^1 = 1_\infty$ if $t > t_1$, $S_t^1 = 1_0$ if $I^{z_1} = 0$, $S_t^1 = 1$ else, and by analogy, $S_t^2 = 2_\infty$ if $t > t_2$, $S_t^2 = 2_0$ if $I^{z_2} = 0$, $S_t^2 = 2$ else.

In what follows, we extend the notation introduced above for partitions of sets. Let us denote $P_t = S_t^1 \dagger S_t^2$ if $t > \tau_{z_1, z_2}^c$ and $P_t = S_t^1 S_t^2$ else. Then, the backward process $(P_t)_{-\infty < t \le t_1}$ is a time in-homogeneous backward Markov chain on the state space $\mathcal{P}_2 = \cup_{a \in \{1_\infty, 1, 1_0\}} \cup_{b \in \{2_\infty, 2, 2_0\}} \{ab, a^\dagger b\}$. At time $t = t_1$, the backward process starts from the initial condition $P_{t_1} = 1 \dagger 2_\infty$ if $t_1 > t_2$, and $P_{t_1} = 1 \dagger 2$ if $t_1 = t_2$ (note that in that case $z_1 \neq z_2$ implies that $i_1 \neq i_2$ so that $P_{t_1} = 12$ is not possible). By Proposition 3.1, the random variable $\min\{\tau_{z_1}^R, \tau_{z_2}^R\} > -\infty$ almost surely, implying that $P_{-\infty} := \lim_{t \to -\infty} P_t$ exists almost surely and $P_{-\infty} \in \{1_0 2_0, 1_0 \dagger 2_0\}$, that is the set of absorbing states. One can check that

$$ (46) \qquad \mathbb{P}_{\theta}(\{z_1 \rightsquigarrow z_2\}) = \mathbb{P}_{\theta}(P_{-\infty} = 1_0 2_0 | P_{t_1} = 1 \dagger 2_{\infty}), $$

in case $t_1 > t_2$, and

$$ (47) \qquad \mathbb{P}_{\theta}(\{z_1 \rightsquigarrow z_2\}) = \mathbb{P}_{\theta}(P_{-\infty} = 1_0 2_0 | P_{t_1} = 1 \dagger 2), $$

in case $t_1 = t_2$.

Now, for $t \le t_1$, let us denote $Q_t(x, y)$ the transition probabilities of the time in-homogeneous backward Markov chain $(P_t)_{-\infty < t \le t_1}$:

$$ Q_t(x, y) = \mathbb{P}_{\theta}(P_{t-1} = y | P_t = x), \quad x, y \in \mathcal{P}_2. $$

Assume for now that $t_1 > t_2$. Since we are only interested in computing the coalescence probability (46), we only need to give the transition probabilities encountered on the path from state $1 \dagger 2_\infty$ to state $1_0 2_0$. For $t > t_2 + 1$, the only relevant transition probability is $Q_t(1 \dagger 2_\infty, 1 \dagger 2_\infty) = 1 - \lambda$. At time $t = t_2 + 1$ there are two relevant transitions: $Q_{t_2+1}(1 \dagger 2_\infty, 1 \dagger 2) = (1 - \lambda)(1 - 1/N)$ and $Q_{t_2+1}(1 \dagger 2_\infty, 12) = (1 - \lambda)/N$.

For $t \le t_2$, there are four relevant transition probabilities:

* $Q_t(1 \dagger 2, 1 \dagger 2) = (1 - \lambda)^2(1 - 1/2)$

* $Q_t(1 \dagger 2, 12) = (1 - \lambda)^2/N$,

* $Q_t(12, 1_0 2_0) = 1 - Q_t(12, 12) = \lambda$.

Figure 3 gives a graphical representation of the relevant transitions just described in case $t_2 < t_1$.

Starting from Equation (46), one can check that

$$ \begin{aligned} \mathbb{P}_{\theta}(\{z_1 \rightsquigarrow z_2\}) = & \mathbb{P}_{\theta}(P_{t_2} = 12 | P_{t_1} = 1 \dagger 2_{\infty}) \mathbb{P}_{\theta}(P_{-\infty} = 1_0 2_0 | P_{t_2} = 12) \\ & + \mathbb{P}_{\theta}(P_{t_2} = 1 \dagger 2 | P_{t_1} = 1 \dagger 2_{\infty}) \mathbb{P}_{\theta}(P_{-\infty} = 1_0 2_0 | P_{t_2} = 1 \dagger 2). \end{aligned} $$

Now, by the dynamics of the Markov chain $(P_t)_{t \le t_1}$ we have that $\mathbb{P}_{\theta}(P_{t_2} = 1 \dagger 2 | P_{t_1} = 1 \dagger 2_\infty) = (1 - \lambda)^{t_1 - t_2}$ and $\mathbb{P}_{\theta}(P_{t_2} = 12 | P_{t_1} = 1 \dagger 2_\infty) = (1 - \lambda)^{t_1 - t_2} / N$. Moreover, $\mathbb{P}_{\theta}(P_{-\infty} = 1_0 2_0 | P_{t_2} = 12) = 1$ by Proposition 3.1. By observing that

$$ \begin{aligned} \mathbb{P}_{\theta}(P_{-\infty} = 1_0 2_0 | P_{t_2} = 1 \dagger 2) &= \frac{(1 - \lambda)^2}{N} \sum_{k=0}^{\infty} (1 - \lambda)^{2k} (1 - 1/N)^k \\ &\le \frac{(1 - \lambda)^2}{N} \frac{1}{1 - (1 - \lambda)^2}, \end{aligned} $$
<figure><img src="image_7.png" /><figcaption>FIG 3. Graph of the relevant transitions of the backward process P used to compute a bound for $IP_{\theta}(\{z_1 \rightsquigarrow z_2\})$ when $t_1 \ge t_2$. The starting node is $1\dagger 2_{\infty}$ if $t_1 > t_2$ and $1\dagger 2$ else. On each edge, the corresponding transition probability is given. The gray vertical line separates two temporal zones: the right one corresponds to times $t > t_2+1$ and the left one corresponds to times $t \le t_2+1$.</figcaption></figure>

and by putting all pieces together, we then obtain that

$$ (48) \quad IP_{\theta}(\{z_1 \rightsquigarrow z_2\}) \le \frac{(1-\lambda)^{|t_1-t_2|}}{N} \left[ 1 + \frac{(1-\lambda)^2}{1-(1-\lambda)^2} \right] = \frac{(1-\lambda)^{|t_1-t_2|}}{N(1-(1-\lambda)^2)}. $$

In the case $t_2 = t_1$, we have

$$ (49) \quad \begin{aligned} IP_{\theta}(\{z_1 \rightsquigarrow z_2\}) &= IP_{\theta}(P_{-\infty} = 1_0 2_0 | P_{t_2} = 1 \dagger 2) \\ &\le \frac{(1-\lambda)^2}{N} \frac{1}{1-(1-\lambda)^2} \le \frac{(1-\lambda)}{N(1-(1-\lambda)^2)}, \end{aligned} $$

so that Item (i) follows from inequalities (48) and (49).

*Proof of Item (ii).* By analogy to the case with two sites, we define $S_t^3 = 3_{\infty}$ if $t > t_3$, $S_t^3 = 3_0$ if $I_t^{z_3} = 0$, $S_t^3 = 3$ else. Let us describe the backward process $(P_t)_{-\infty<t\le t_1}$ in that case. First note that for all $t \le t_1$, $P_t$ is a partition of $\{S_t^1, S_t^2, S_t^3\}$ and in particular $P_s \in \mathcal{P}_3 = \bigcup_{a \in \{1_{\infty}, 1, 1_0\}} \bigcup_{b \in \{2_{\infty}, 2, 2_0\}} \bigcup_{c \in \{3_{\infty}, 3, 3_0\}} \mathcal{P}(\{a, b, c\})$. The choice of the partition is induced by the equivalence relation defined by

$$ S_t^k \leftrightarrow_t S_t^\ell \text{ if and only if } t \le \tau_{z_k, z_\ell}^c. $$

This equivalence relation naturally induces a partition of $\{S_t^1, S_t^2, S_t^3\}$ and we define $P_t$ as this partition. Similarly to the proof of Item (i), we have

$$ (50) \quad IP_{\theta}(\{z_1 \rightsquigarrow z_2 \rightsquigarrow z_3\}) = IP_{\theta}(P_{-\infty} = 1_0 2_0 3_0 | P_{t_1} = 1 \dagger 2_{\infty} \dagger 3_{\infty}), $$

in case $t_1 > t_2 \ge t_3$ for instance, where $P_{-\infty} = \lim_{t \to -\infty} P_t$ exists almost surely and $P_{-\infty} \in \mathcal{P}(\{1_0, 2_0, 3_0\})$. Assume for now that $t_1 > t_2 > t_3$ which is the most general case. The other cases can be treated similarly.

For $t \le t_1$, let us denote the transition probabilities of the chain

$$ Q_t(x, y) = IP_{\theta}(P_{t-1} = y | P_t = x), x, y \in \mathcal{P}_3. $$
Since we are only interested in computing the coalescence probability (50), we only need to give the transition probabilities encountered on the path from state $1\dagger 2_{\infty}\dagger 3_{\infty}$ to state $1_02_03_0$. For $t > t_2 + 1$, the only relevant transition probability is $Q_t(1\dagger 2_{\infty}\dagger 3_{\infty}, 1\dagger 2_{\infty}\dagger 3_{\infty}) = 1 - \lambda$. At time $t = t_2 + 1$, there are two relevant transition probabilities: $Q_{t_2+1}(1\dagger 2_{\infty}\dagger 3_{\infty}, 1\dagger 2\dagger 3_{\infty}) = (1 - \lambda)(1 - 1/N)$ and $Q_{t_2+1}(1\dagger 2_{\infty}\dagger 3_{\infty}, 12\dagger 3_{\infty}) = (1 - \lambda)/(N)$.

For $t_3 + 1 < t \le t_2$ the three relevant transition probabilities are:

* $Q_t(1\dagger 2\dagger 3_{\infty}, 1\dagger 2\dagger 3_{\infty}) = (1 - \lambda)^2(1 - 1/N) \le (1 - \lambda)^2$,

* $Q_t(1\dagger 2\dagger 3_{\infty}, 12\dagger 3_{\infty}) = (1 - \lambda)^2/N \le (1 - \lambda)N^{-1}$,

* $Q_t(12\dagger 3_{\infty}, 12\dagger 3_{\infty}) = (1 - \lambda)$.

At time $t = t_3 + 1$, the seven relevant transition probabilities are:

* $Q_t(1\dagger 2\dagger 3_{\infty}, 1\dagger 2\dagger 3) = (1 - \lambda)^2(1 - 1/N)(1 - 2/N) \le (1 - \lambda)^2$,

* $Q_t(1\dagger 2\dagger 3_{\infty}, 1\dagger 23) = Q_t(1\dagger 2\dagger 3_{\infty}, 13\dagger 2) = Q_t(1\dagger 2\dagger 3_{\infty}, 12\dagger 3) \le (1 - \lambda)^2 N^{-1}$,

* $Q_t(1\dagger 2\dagger 3_{\infty}, 123) = (1 - \lambda)^2/N^2$.

* $Q_t(12\dagger 3_{\infty}, 12\dagger 3) = (1 - \lambda)(1 - N) \le (1 - \lambda)$,

* $Q_t(12\dagger 3_{\infty}, 123) = (1 - \lambda)/N$.

Finally, for $t \le t_3$, the thirteen relevant transition probabilities are:

* $Q_t(1\dagger 2\dagger 3, 1\dagger 2\dagger 3) \le (1 - \lambda)^3$,

* $Q_t(1\dagger 2\dagger 3, 12\dagger 3) = Q_t(1\dagger 2\dagger 3, 1\dagger 23) = Q_t(1\dagger 2\dagger 3, 2\dagger 13) \le (1 - \lambda)^3 N^{-1}$,

* $Q_t(1\dagger 2\dagger 3, 123) \le (1 - \lambda)^3 N^{-2}$,

* $Q_t(12\dagger 3, 12\dagger 3) = Q_t(1\dagger 23, 1\dagger 23) = Q_t(2\dagger 13, 2\dagger 13) \le (1 - \lambda)^2$,

* $Q_t(12\dagger 3, 123) = Q_t(1\dagger 23, 123) = Q_t(2\dagger 13, 123) \le (1 - \lambda)^2 N^{-1}$,

* $Q_t(123, 1_02_03_0) = 1 - Q_t(123, 123) = \lambda$.

Looking at the Figures 4 and 5, the Markov property of $(P_t)_{t \le t_1}$ permits to compute the probabilities of ending in state $1_02_03_0$ from any of the other states of the graph. These two figures highlight three time zones denoted I, II and III, which are depicted from left to right corresponding to the forward evolution of time (hence the arrows go from right to left since they correspond to a backward dynamics).

In the following, we decompose the computation of the probability in Equation (50) with respect to the three zones. To simplify the notation, let us denote for all $t \in \mathbb{Z}$ and $\Theta \in \mathcal{P}_3$, $p(t, \Theta) = \mathbb{P}_{\Theta}(P_{-\infty} = 1_02_03_0 | P_t = \Theta)$.

For zone I, we have to consider initial conditions starting at time $t_3$ (see Figure 5). That is,

$$p(t_3, 1\dagger 23) = p(t_3, 12\dagger 3) = p(t_3, 13\dagger 2) \le K(1 - \lambda)^2 N^{-1},$$

and

$$p(t_3, 1\dagger 2\dagger 3) \le K(1 - \lambda)^3 [3p(t_3, 1\dagger 23)N^{-1} + N^{-2}] \le K(1 - \lambda)^3 N^{-2}.$$

For zone II, we have for all $t_3 < t \le t_2$,

$$p(t, 12\dagger 3_{\infty}) = (1 - \lambda)^{t-t_3} [p(t_3, 12\dagger 3) + N^{-1}p(t_3, 123)] \le K(1 - \lambda)^{t-t_3} N^{-1},$$

since we have to wait in the middle part during $t - t_3 - 1$ steps and then take the edge from $12\dagger 3_{\infty}$ to $12\dagger 3$ or the edge from $12\dagger 3_{\infty}$ to $123$. Moreover, we have

$$p(t_2, 1\dagger 2\dagger 3_{\infty}) \le K \left\{ \begin{aligned} & \sum_{k=1}^{t_2-t_3-1} (1-\lambda)^{2k} p(t_2-k, 12\dagger 3_{\infty}) N^{-1} \\ & + (1-\lambda)^{2(t_2-t_3)} [3p(t_3, 1\dagger 23) N^{-1} + p(t_3, 1\dagger 2\dagger 3)] \end{aligned} \right\}$$
<figure><img src="image_8.png" /><figcaption>FIG 4. Graph of the relevant transitions of the backward process P (for times $t > t_3 + 1$) used to compute a bound for $\mathbb{P}_{\theta}(\{z_1 \rightsquigarrow z_2 \rightsquigarrow z_3\})$ when $t_1 \ge t_2 \ge t_3$. The starting node is $1\dagger 2_\infty \dagger 2_\infty$ if $t_1 > t_2 > 3$, $1\dagger 2\dagger 3_\infty$ if $t_1 = t_2 > t_3$ and $1\dagger 2\dagger 3$ (which appears in Figure 5) else. On each edge, an upper bound of the corresponding transition probability is given. The gray vertical line separates two temporal zones: zone II corresponds to times $t_3 + 1 < t \le t_2 + 1$ and zone III corresponds to times $t > t_2 + 1$.</figcaption></figure>

$$ \begin{aligned} & \le K \left\{ (1-\lambda)^{t_2-t_3} \sum_{k=1}^{t_2-t_3-1} (1-\lambda)^k N^{-2} + 4(1-\lambda)^{2(t_2-t_3)} (1-\lambda)^2 N^{-2} \right\} \\ & \le K(1-\lambda)^{t_2-t_3} N^{-2}. \end{aligned} $$

Finally, using zone III, we have

$$ \begin{aligned} \mathbb{P}_{\theta}(\{z_1 \rightsquigarrow z_2 \rightsquigarrow z_3\}) &= p(t_1, 1\dagger 2_{\infty} \dagger 3_{\infty}) \\ &= (1-\lambda)^{t_1-t_2} [p(t_2, 1\dagger 2\dagger 3_{\infty})(1-1/2N) + p(t_2, 12\dagger 3_{\infty})(2N)^{-1}], \end{aligned} $$

which ends the proof.

**Proof of Item (iii).** The backward process $(P_t)_{-\infty<t\le t_1}$ is naturally extended to the case of four sites and we use the decomposition

$$ \{z_1 \rightsquigarrow z_2\} \cap \{z_3 \rightsquigarrow z_4\} = \{z_1 \rightsquigarrow z_2 \rightsquigarrow z_3 \rightsquigarrow z_4\} \cup \{z_1 \rightsquigarrow z_2 \rightsquigarrow z_3 \rightsquigarrow z_4\}. $$

Assume for now that $t_1 > t_2 \ge t_3 \ge t_4$ so that $P_{t_1} = 1\dagger 2_\infty \dagger 3_\infty \dagger 4_\infty$. It suffices to control

$$ \begin{aligned} \mathbb{P}_{\theta}(P_{-\infty} = 1_0 2_0 3_0 4_0 | P_{t_1} = 1\dagger 2_{\infty} \dagger 3_{\infty} \dagger 4_{\infty}) \\ + \mathbb{P}_{\theta}(P_{-\infty} = 1_0 2_0 \dagger 3_0 4_0 | P_{t_1} = 1\dagger 2_{\infty} \dagger 3_{\infty} \dagger 4_{\infty}). \end{aligned} $$

The list of the relevant transitions and the corresponding graphs are not written here due to their huge size. However the computations can be made in the same manner. For instance, one can prove that

$$ \mathbb{P}_{\theta}(P_{-\infty} = 1_0 2_0 3_0 4_0 | P_{t_1} = 1\dagger 2_{\infty} \dagger 3_{\infty} \dagger 4_{\infty}) \le K(1-\lambda)^{t_1-t_4} N^{-3}, $$
<figure><img src="image_9.png" /><figcaption>FIG 5. Graph of the relevant transitions of the backward process P (for times $t \le t_3 + 1$) used to compute a bound for $IP_\theta(\{z_1 \rightsquigarrow z_2 \rightsquigarrow z_3\})$ when $t_1 \ge t_2 \ge t_3$. The starting node is $1^\dagger 2_\infty \dagger 2_\infty$ (which appears in Figure 4) if $t_1 > t_2 > 3$, $1^\dagger 2 \dagger 3_\infty$ if $t_1 = t_2 > t_3$ and $1^\dagger 2 \dagger 3$ else. On some of the edges, an upper bound of the corresponding transition probability is given (the nodes $1^\dagger 2$, $13^\dagger 2$ and $12^\dagger 3$ play almost the same role that is why we omit some the transition probabilities). The gray vertical line separates two temporal zones: zone I corresponds to times $t \le t_3 + 1$ and zone II corresponds to times $t_3 + 1 < t \le t_2 + 1$.</figcaption></figure>

and

$$IP_\theta(P_{-\infty} = 1_0 2_0 \dagger 3_0 4_0 | P_{t_1} = 1^\dagger 2_\infty \dagger 3_\infty \dagger 4_\infty) \le K(1-\lambda)^{(t_1-t_2)+(t_3-t_4)} N^{-2},$$

which in turn give the desired result. Finally, the same upper bound holds if $t_1 > t_2$ is not satisfied by replacing the initial condition (for instance $t_1 = t_2 > t_3 \ge t_4$ corresponds to the initial condition $P_{t_1} = 1^\dagger 2^\dagger 3_\infty \dagger 4_\infty$ if $z_1 \neq z_2$).

*Proof of Item (iv).* We will only prove the case $t_1 \wedge t_2 > t_3 \vee t_4$. The case in which $t_1 \wedge t_2 = t_3 \vee t_4$ follows along the same lines. In the sequel, as before, we assume without loss of generality that $t_1 \ge t_2 > t_3 \ge t_4$ (the case where $t_1 < t_2$ can be treated in the same way). In the following, we use the backward process $(P_t)_{-\infty<t\le t_1}$ with four sites like in the proof of Item (iii). Then, since $t_3 + 1 \le t_2$,

$$ \begin{aligned} & \{z_1 \rightsquigarrow z_3 \rightsquigarrow z_4 \not\rightsquigarrow z_2\} \cap \{\tau_{z_2}^R < t_3\} \subset \{z_1 \rightsquigarrow z_3 \rightsquigarrow z_4 \not\rightsquigarrow z_2\} \cap \{\tau_{z_2}^R < t_3 + 1\} \\ & = \{P_{-\infty} = 1_0 3_0 4_0 \dagger 2_0\} \cap \{P_{t_3+1} = 1^\dagger 2^\dagger 3_\infty \dagger 4_\infty\}. \end{aligned} $$
Assume for now that $t_1 > t_2$. It follows from the same arguments as above that

$$ \mathbb{P}_{\theta}(P_{-\infty} = 1_03_04_0\dagger 2_0 | P_{t_3+1} = 1\dagger 2\dagger 3_{\infty}\dagger 4_{\infty}) \le K(1-\lambda)^{t_3+1-t_4}N^{-2}, $$

and

$$ \mathbb{P}_{\theta}(P_{t_3+1} = 1\dagger 2\dagger 3_{\infty}\dagger 4_{\infty} | P_{t_1} = 1\dagger 2_{\infty}\dagger 3_{\infty}\dagger 4_{\infty}) \le (1-\lambda)^{t_1-(t_3+1)+t_2-(t_3+1)}. $$

which in turn give the desired result.

Finally, if $t_1 = t_2$, the same upper-bound holds with the initial condition $P_{t_1} = 1\dagger 2_{\infty}\dagger 3_{\infty}\dagger 4_{\infty}$ replaced by $P_{t_1} = 1\dagger 2\dagger 3_{\infty}\dagger 4_{\infty}$. $\square$

## APPENDIX D: PROOF OF LEMMA 3.8

The statements in both Items 1 and 2 are trivial, so that we only need to prove Items 3 through 7. Items 3, 4 and 5 are proved in a similar way: Lemma B.1 is used to relate independence with the high probability event of no coalescence, and some easy computations are needed to bound the covariance. Items 6 and 7 use the more involved Lemma B.2.

### PROOF OF ITEM 3.
In this case, $B = \text{Cov}_{\theta}[Y_{z_1}^2, Y_{z_3}^2]$. Let $\tilde{X}_{z_1}$ be a random variable distributed as $X_{z_1}$, independent of $X_{z_3}$ and such that $\{X_{z_1} \ne \tilde{X}_{z_1}\} \subset \{z_1 \rightsquigarrow z_3\}$. The existence of such random variables is ensured by Lemma B.1 with $k=2$. In particular, we have that $\mathbb{E}_{\theta}[\tilde{X}_{z_1}] = \mathbb{E}_{\theta}[X_{z_1}]$. Denote $\tilde{Y}_{z_1} = \tilde{X}_{z_1} - \mathbb{E}_{\theta}[\tilde{X}_{z_1}]$ and observe that $\{\tilde{Y}_{z_1} \ne Y_{z_1}\} = \{\tilde{X}_{z_1} \ne X_{z_1}\}$. By using the properties satisfied by random variables $\tilde{X}_{z_1}$, one can check that $\tilde{Y}_{z_1}$ is independent of $Y_{z_3}$, so that

$$ B = \text{Cov}_{\theta}[Y_{z_1}^2, Y_{z_3}^2] = \mathbb{E}_{\theta}[(Y_{z_1}^2 - \tilde{Y}_{z_1}^2)Y_{z_3}^2]. $$

Since $\{\tilde{Y}_{z_1} \ne Y_{z_1}\} = \{\tilde{X}_{z_1} \ne X_{z_1}\} \subset \{z_1 \rightsquigarrow z_3\}$ and $\max\{|Y_{z_1}|, |\tilde{Y}_{z_1}|, |Y_{z_3}|\} \le 1$ almost surely, it follows that

$$ |B| \le \mathbb{P}_{\theta}(z_1 \rightsquigarrow z_3), $$

so that the result follows Item (i) of Proposition 3.5. $\square$

<figure><img src="image_10.png" /></figure>

### PROOF OF ITEM 4.
Suppose that $z_1 = z_2$. Then $B = \text{Cov}_{\theta}[Y_{z_1}^2, Y_{z_3}Y_{z_4}]$. By Lemma B.1 with $k=3$, there exist random variables $\tilde{Y}_{z_1}, \tilde{Y}_{z_3}$ and $\tilde{Y}_{z_4}$ defined in such a way that the following properties hold for each $i \in \{1, 3, 4\}$: 1) $\tilde{Y}_{z_i}$ has the same law as $Y_{z_i}$; 2) $\tilde{Y}_{z_i}$ is independent of $\tilde{Y}_{z_j}$ and $Y_{z_j}$ for all $j \in \{1, 3, 4\} \setminus \{i\}$; 3) $\{\tilde{Y}_{z_i} \ne Y_{z_i}\} \subset \cup_{j \in \{1, 3, 4\} \setminus \{i\}} \{z_i \rightsquigarrow z_j\}$. Using first that $Y_{z_1}^2$ has same law as $\tilde{Y}_{z_1}^2$ and then that $\mathbb{E}_{\theta}[\tilde{Y}_{z_i}] = \mathbb{E}_{\theta}[Y_{z_i}] = 0$, we can deduce from properties 1 and 2 above that

$$ \begin{aligned} B &= \mathbb{E}_{\theta}[(Y_{z_1}^2 - \tilde{Y}_{z_1}^2)Y_{z_3}Y_{z_4}] \\ &= \mathbb{E}_{\theta}[(Y_{z_1}^2 - \tilde{Y}_{z_1}^2)(Y_{z_3} - \tilde{Y}_{z_3})(Y_{z_4} - \tilde{Y}_{z_4})]. \end{aligned} $$

By combining property 3, the above identity and the fact that $\rightsquigarrow$ is an equivalence relation, we can deduce that

$$ |B| \le \mathbb{E}_{\theta}[1_{\{z_1 \rightsquigarrow z_3 \rightsquigarrow z_4\}}] = \mathbb{P}_{\theta}(\{z_1 \rightsquigarrow z_2 \rightsquigarrow z_3\}), $$

so that the result follows from Item (ii) of Proposition 3.5. $\square$

<figure><img src="image_11.png" /></figure>

### PROOF OF ITEM 5.
Let us consider the case in which $z_1 = z_3$. Let $z_1, z_2, z_4$ be three different sites. The goal is to bound the covariance

$$ (51) \quad B = \text{Cov}_{\theta}[Y_{z_1}Y_{z_2}, Y_{z_1}Y_{z_4}] = \mathbb{E}_{\theta}[Y_{z_1}^2 Y_{z_2} Y_{z_4}] - \text{Cov}_{\theta}[Y_{z_1}, Y_{z_2}] \text{Cov}_{\theta}[Y_{z_1}, Y_{z_4}]. $$
Arguing exactly like in the proof of Item 3, one can show that $\mathrm{Cov}_{\theta}[Y_{z_1}, Y_{z_2}] \le K(1 - \lambda)^{|t_1 - t_2|} N^{-1}$ and $\mathrm{Cov}_{\theta}[Y_{z_1}, Y_{z_4}] \le K(1 - \lambda)^{|t_1 - t_4|} N^{-1}$, so that

$$ (52) \qquad \mathrm{Cov}_{\theta}[Y_{z_1}, Y_{z_2}] \mathrm{Cov}_{\theta}[Y_{z_1}, Y_{z_4}] \le K(1 - \lambda)^{|t_1 - t_2| + |t_1 - t_4|} N^{-2}. $$

We will now deal with the term $\mathbb{E}_{\theta}[Y_{z_1}^2 Y_{z_2} Y_{z_4}]$. For that end, let $\tilde{X}_{z_1}$, $\tilde{X}_{z_2}$ and $\tilde{X}_{z_4}$ be the random variables defined in Lemma B.1 with $k=3$, and denote $\tilde{Y}_{z_i} = \tilde{X}_{z_i} - \mathbb{E}_{\theta}(X_{z_i})$, $i \in \{1, 2, 4\}$ the corresponding centered versions. By using the independence properties of the random variables $\tilde{Y}_{z_i}$ (and the fact that they are all centered), one can check that

$$ \mathbb{E}_{\theta}[Y_{z_1}^2 Y_{z_2} Y_{z_4}] = \mathbb{E}_{\theta}[(Y_{z_1}^2 - \tilde{Y}_{z_1}^2)(Y_{z_2} - \tilde{Y}_{z_2})(Y_{z_4} - \tilde{Y}_{z_4})] + \mathbb{E}_{\theta}[\tilde{Y}_{z_1}^2] \mathrm{Cov}_{\theta}[Y_{z_2}, Y_{z_4}]. $$

Arguing exactly like in the proof of Item 4, one can deduce that

$$ (53) \qquad \mathbb{E}_{\theta}[(Y_{z_1}^2 - \tilde{Y}_{z_1}^2)(Y_{z_2} - \tilde{Y}_{z_2})(Y_{z_4} - \tilde{Y}_{z_4})] \le K(1 - \lambda)^{|s_3 - s_2| + |s_2 - s_1|} N^{-2}, $$

where $s_1 \le s_2 \le s_3$ is an ordering of $\{t_1, t_2, t_4\}$. Finally, since $\mathrm{Cov}_{\theta}[Y_{z_2}, Y_{z_4}] \le K(1 - \lambda)^{|t_2 - t_4|} N^{-1}$, the random variables $\tilde{Y}_{z_k}$ are bounded by 1 almost surely and $|t_2 - t_4| \le \min\{|t_2 - t_1| + |t_1 - t_4|, |s_3 - s_2| + |s_2 - s_1|\}$, it follows from (51), (52) and (53) that

$$ |B| \le K(1 - \lambda)^{|t_2 - t_4|} N^{-1}, $$

concluding the proof.

<figure><img src="image_12.png" /></figure>

PROOF OF ITEM 6. Let $z_1, \dots, z_4$ be four different points in $\mathcal{Z}$ and denote by $\hat{X}_{z_1}, \hat{X}_{z_2}, \tilde{X}_{z_1}, \dots, \tilde{X}_{z_4}$ the random variables defined in Lemma B.2. Moreover, denote by $\tilde{Y}_{z_k} = \hat{X}_{z_k} - \mathbb{E}_{\theta}(\tilde{X}_{z_k})$, $1 \le k \le 4$ and $\hat{Y}_{z_k} = \hat{X}_{z_k} - \mathbb{E}_{\theta}(\hat{X}_{z_k})$, $1 \le k \le 2$, their centered versions. First of all, observe that by using Items (i), (ii) and (iii) of Lemma B.2, we can rewrite B as

$$ B = \mathbb{E}_{\theta}\left[(Y_{z_1}Y_{z_2} - \hat{Y}_{z_1}\hat{Y}_{z_2})Y_{z_3}Y_{z_4}1_{\{z_1 \rightsquigarrow z_3\} \cup \{z_1 \rightsquigarrow z_4\} \cup \{z_2 \rightsquigarrow z_3\} \cup \{z_2 \rightsquigarrow z_4\}}\right]. $$

In the remaining of the proof, we adopt the following notation. For $i, j, k, \ell \in \{1, \dots, 4\}$, we write $\{ijk\ell\}$, $\{ijk|\ell\}$, $\{ij|k\ell\}$ and $\{ij|k|\ell\}$ to denote, respectively, the events $\{z_i \rightsquigarrow z_j \rightsquigarrow z_k\}$, $\{z_i \rightsquigarrow z_j \rightsquigarrow z_k \nparallel z_\ell\}$, $\{z_i \rightsquigarrow z_j \nparallel z_k \nparallel z_\ell\}$ and $\{z_k \nparallel z_i \rightsquigarrow z_j \nparallel z_\ell\}$. With this notation, observe that we can write

$$ 1_{\{z_1 \rightsquigarrow z_3\} \cup \{z_1 \rightsquigarrow z_4\} \cup \{z_2 \rightsquigarrow z_3\} \cup \{z_2 \rightsquigarrow z_4\}} = 1_{\{1234\}} + 1_{\{123|4\}} + 1_{\{124|3\}} + 1_{\{134|2\}} + 1_{\{234|1\}} \\ + 1_{\{13|24\}} + 1_{\{14|23\}} + 1_{\{13|2|4\}} + 1_{\{14|2|3\}} + 1_{\{23|1|4\}} + 1_{\{24|1|4\}}. $$

Let us denote

$$ M = \max\{\mathbb{P}_{\theta}(1234), \mathbb{P}_{\theta}(13|24), \mathbb{P}_{\theta}(14|23), \\ \mathbb{P}_{\theta}(\{134|2\} \cap \{\tau_{z_2}^R < t_3 \lor t_4\}), \mathbb{P}_{\theta}(\{234|1\} \cap \{\tau_{z_1}^R < t_3 \lor t_4\})\}. $$

Step 1. Since the random variables $Y_{z_i}$'s and $\hat{Y}_{z_i}$'s are bounded by 1 almost surely, we clearly have

$$ \left| \mathbb{E}_{\theta} \left[ \left( Y_{z_1} Y_{z_2} - \hat{Y}_{z_1} \hat{Y}_{z_2} \right) Y_{z_3} Y_{z_4} \left( 1_{\{1234\}} + 1_{\{13|24\}} + 1_{\{14|23\}} \right) \right] \right| \le 6M. $$
Step 2. Let E denote one of the following events $\{123|4\}$ or $\{124|3\}$. Then, we prove the following inequality

$$ (54) \qquad \left| \mathbb{E}_{\theta} \left[ \left( Y_{z_1} Y_{z_2} - \hat{Y}_{z_1} \hat{Y}_{z_2} \right) Y_{z_3} Y_{z_4} 1_E \right] \right| \le 2M. $$

To do so, let us assume that $E = \{123|4\}$. The other case is treated similarly. By Item (vi) of Lemma B.2, we have that $Y_{z_4} = \tilde{Y}_{z_4}$ on E, so that

$$ \mathbb{E}_{\theta} \left[ \left( Y_{z_1} Y_{z_2} - \hat{Y}_{z_1} \hat{Y}_{z_2} \right) Y_{z_3} Y_{z_4} 1_E \right] = \mathbb{E}_{\theta} \left[ \left( Y_{z_1} Y_{z_2} - \hat{Y}_{z_1} \hat{Y}_{z_2} \right) Y_{z_3} \tilde{Y}_{z_4} 1_E \right]. $$

By using that $1_{\{123|4\}} + 1_{\{1234\}} = 1_{\{123\}}$ is $\sigma(I^{z_1}, I^{z_2}, I^{z_3})$-measurable and the fact that $\tilde{Y}_{z_4}$ is a centered random variable which is independent of $(Y_{z_1} Y_{z_2} - \hat{Y}_{z_1} \hat{Y}_{z_2}) Y_{z_3} 1_{\{123\}}$ (thanks to Items (ii) and (v) of Lemma B.2), it follows that

$$ \mathbb{E}_{\theta} \left[ \left( Y_{z_1} Y_{z_2} - \hat{Y}_{z_1} \hat{Y}_{z_2} \right) Y_{z_3} \tilde{Y}_{z_4} 1_E \right] = - \mathbb{E}_{\theta} \left[ \left( Y_{z_1} Y_{z_2} - \hat{Y}_{z_1} \hat{Y}_{z_2} \right) Y_{z_3} \tilde{Y}_{z_4} 1_{\{1234\}} \right]. $$

Finally, we conclude like in Step 1.

Step 3. Let E denote one of the following events $\{134|2\}$ or $\{234|1\}$. Then, we prove the following inequality

$$ (55) \qquad \left| \mathbb{E}_{\theta} \left[ \left( Y_{z_1} Y_{z_2} - \hat{Y}_{z_1} \hat{Y}_{z_2} \right) Y_{z_3} Y_{z_4} 1_E \right] \right| \le 5M. $$

To do so, let us assume that $E = \{134|2\}$. The other case is treated similarly. In the following, let us denote $E_1 = E \cap \{\tau_{z_2}^R \ge t_3 \lor t_4\}$ and $E_2 = E \cap \{\tau_{z_2}^R < t_3 \lor t_4\}$. On the one hand, we use the fact that both

$$ \mathbb{E}_{\theta} \left[ \left( Y_{z_1} Y_{z_2} - \hat{Y}_{z_1} \hat{Y}_{z_2} \right) Y_{z_3} Y_{z_4} 1_{E_2} \right] \text{ and } \mathbb{E}_{\theta} \left[ \left( Y_{z_1} - \hat{Y}_{z_1} \right) \tilde{Y}_{z_2} Y_{z_3} Y_{z_4} 1_{E_2} \right] $$

are upper bounded by $2IP_{\theta}(E_2)$. On the other hand, we know by Items (iii) and (vi) of Lemma B.2, that $\tilde{Y}_{z_2} = Y_{z_2} = \hat{Y}_{z_2}$ on $E_1$. Combining those two properties, we can get

$$ \left| \mathbb{E}_{\theta} \left[ \left( Y_{z_1} Y_{z_2} - \hat{Y}_{z_1} \hat{Y}_{z_2} \right) Y_{z_3} Y_{z_4} 1_E \right] \right| \le \left| \mathbb{E}_{\theta} \left[ \left( Y_{z_1} - \hat{Y}_{z_1} \right) \tilde{Y}_{z_2} Y_{z_3} Y_{z_4} 1_E \right] \right| + 4IP_{\theta}(E_2). $$

By using that $1_{\{134|2\}} + 1_{\{1234\}} = 1_{\{134\}}$ is $\sigma(I^{z_1}, I^{z_3}, I^{z_4})$-measurable and the fact that $\tilde{Y}_2$ is a centered random variable which is independent of $(Y_{z_1} - \hat{Y}_{z_1}) Y_{z_3} Y_{z_4} 1_{\{134\}}$ (thanks to Item (v) of Lemma B.2), it follows that

$$ \mathbb{E}_{\theta} \left[ \left( Y_{z_1} - \hat{Y}_{z_1} \right) \tilde{Y}_{z_2} Y_{z_3} Y_{z_4} 1_{\{134|2\}} \right] = - \mathbb{E}_{\theta} \left[ \left( Y_{z_1} - \hat{Y}_{z_1} \right) \tilde{Y}_{z_2} Y_{z_3} Y_{z_4} 1_{\{1234\}} \right]. $$

Finally, we conclude like in Step 1.

Step 4. Let E denote one of the following events: $\{13|2|4\}$, $\{14|2|3\}$, $\{23|1|4\}$ or $\{24|1|3\}$. Then, we prove the following inequality

$$ \left| \mathbb{E}_{\theta} \left[ \left( Y_{z_1} Y_{z_2} - \hat{Y}_{z_1} \hat{Y}_{z_2} \right) Y_{z_3} Y_{z_4} 1_E \right] \right| \le 7M. $$

To do so, let us assume that $E = \{13|2|4\}$. The other cases are treated similarly. First, we combine the fact that $Y_{z_4} = \tilde{Y}_{z_4}$ on E (this holds by Item (vi) of Lemma B.2) together with the decomposition $1_{\{13|2|4\}} = 1_{13|2} - 1_{13|24} - 1_{134|2}$ to deduce that

$$ \begin{aligned} \mathbb{E}_{\theta} \left[ \left( Y_{z_1} Y_{z_2} - \hat{Y}_{z_1} \hat{Y}_{z_2} \right) Y_{z_3} Y_{z_4} 1_{\{13|2|4\}} \right] = & \mathbb{E}_{\theta} \left[ \left( Y_{z_1} Y_{z_2} - \hat{Y}_{z_1} \hat{Y}_{z_2} \right) Y_{z_3} \tilde{Y}_{z_4} 1_{\{13|2\}} \right] \\ & - \mathbb{E}_{\theta} \left[ \left( Y_{z_1} Y_{z_2} - \hat{Y}_{z_1} \hat{Y}_{z_2} \right) Y_{z_3} \tilde{Y}_{z_4} 1_{\{13|24\}} \right] - \mathbb{E}_{\theta} \left[ \left( Y_{z_1} Y_{z_2} - \hat{Y}_{z_1} \hat{Y}_{z_2} \right) Y_{z_3} \tilde{Y}_{z_4} 1_{\{134|2\}} \right]. \end{aligned} $$
Now, since $\tilde{Y}_{z_4}$ is centered and independent of $Y_{z_1}, \hat{Y}_{z_1}, Y_{z_2}, \hat{Y}_{z_2}, Y_{z_3}$ and $1_{13|2}$ (because this random variable is $\sigma(I^{z_1}, I^{z_2}, I^{z_3})$-measurable), we have that

$$ \mathbb{E}_{\theta} \left[ \left( Y_{z_1} Y_{z_2} - \hat{Y}_{z_1} \hat{Y}_{z_2} \right) Y_{z_3} \tilde{Y}_{z_4} 1_{\{13|2\}} \right] = 0, $$

which implies that

$$ \begin{aligned} \mathbb{E}_{\theta} \left[ \left( Y_{z_1} Y_{z_2} - \hat{Y}_{z_1} \hat{Y}_{z_2} \right) Y_{z_3} Y_{z_4} 1_{\{13|2|4\}} \right] = & - \mathbb{E}_{\theta} \left[ \left( Y_{z_1} Y_{z_2} - \hat{Y}_{z_1} \hat{Y}_{z_2} \right) Y_{z_3} \tilde{Y}_{z_4} 1_{\{13|24\}} \right] \\ & - \mathbb{E}_{\theta} \left[ \left( Y_{z_1} Y_{z_2} - \hat{Y}_{z_1} \hat{Y}_{z_2} \right) Y_{z_3} \tilde{Y}_{z_4} 1_{\{134|2\}} \right]. \end{aligned} $$

The first term on the right-hand side of the above equality is in absolute value at most $2M$. The second one can be dealt with proceeding similarly as Step 3.

Combining Steps 1 to 4, we get $|B| \le 48M$, and Items (iii) and (iv) of Proposition 3.5 give the upper-bound $M \le K(1 - \lambda)^{t_1 \vee t_2 - t_3 \wedge t_4} N^{-2}$, which allows to conclude. $\square$

PROOF OF ITEM 7. By symmetry, we only need to consider the case $t_1 \le t_2$ and $t_3 \le t_4$. The case $t_1 \ge t_3$ or $t_2 \le t_3$ has already been treated in Item 6. So suppose that either $t_1 \le t_3 \le t_2 \le t_4$ or that $t_3 \le t_1 \le t_4 \le t_2$. Since $\mathbb{E}_{\theta}[Y_{z_i}] = 0$, we can write

$$ (56) \quad B = \mathrm{Cov}_{\theta}[Y_{z_1} Y_{z_3}, Y_{z_2} Y_{z_4}] + \mathrm{Cov}_{\theta}[Y_{z_1}, Y_{z_3}] \mathrm{Cov}_{\theta}[Y_{z_2}, Y_{z_4}] \\ \qquad - \mathrm{Cov}_{\theta}[Y_{z_1}, Y_{z_2}] \mathrm{Cov}_{\theta}[Y_{z_3}, Y_{z_4}]. $$

By Lemma 3.7, it follows that

$$ (57) \quad \begin{aligned} & |\mathrm{Cov}_{\theta}[Y_{z_1}, Y_{z_3}] \mathrm{Cov}_{\theta}[Y_{z_2}, Y_{z_4}]| + |\mathrm{Cov}_{\theta}[Y_{z_1}, Y_{z_2}] \mathrm{Cov}_{\theta}[Y_{z_3}, Y_{z_4}]| \\ & \le KN^{-2} \left[ (1-\lambda)^{|t_1-t_3|+|t_2-t_4|} + (1-\lambda)^{|t_1-t_2|+|t_3-t_4|} \right]. \end{aligned} $$

Also, from Item 6 above, we have that

$$ (58) \quad \mathrm{Cov}_{\theta}[Y_{z_1} Y_{z_3}, Y_{z_2} Y_{z_4}] \le K(1-\lambda)^{\max t_i - \min t_i} N^{-2}. $$

Combining (56), (57) and (58), we obtain the assertion. $\square$

APPENDIX E: TEMPORAL CONVERGENCES

**E.1. Proof of Equation (32).** This one is simple. Since we start from the stationary distribution, we have

$$ \mathbb{E}_{\theta}[\hat{m}] = T^{-1}N^{-1} \sum_{t=1}^{T} \sum_{i=1}^{N} \mathbb{E}_{\theta}[X_{i,t}] = N^{-1} \sum_{i=1}^{N} m_i^N = \overline{m^N}. $$

Hence, the left-hand side of Equation (32) is equal to $\mathrm{Var}_{\theta}(\hat{m})$.

Then, note that $\hat{m} = (TN)^{-1} \sum_{t=1}^{T} \sum_{i=1}^{N} X_{i,t}$ so that

$$ \begin{aligned} \mathrm{Var}_{\theta}(\hat{m}) &\le (TN)^{-2} \sum_{t_1,t_2=1}^{T} \sum_{i_1,i_2=1}^{N} |\mathrm{Cov}_{\theta}(X_{i_1,t_1}, X_{i_2,t_2})| \\ &\le (TN)^{-1} \left[ \frac{2K}{\lambda} + K(1-\lambda) + 1 \right], \end{aligned} $$

where the last inequality comes from Lemma 3.7.
**E.2. Proof of Equation (34).** The proof relies on a decomposition of the error $|\hat{w} - w_{\infty}^N|$ (see Inequality (59)). Each term of the decomposition is then controlled by one of the three Propositions E.1, E.7 or E.10. Their proofs need intermediate results, which are provided as lemmas. They are linked with martingale properties of the process (Lemmas E.2, E.3 and E.5) or with the vanishing covariances (Lemma E.8).

Recall that

$$ \hat{w} = 2W_{2\Delta} - W_{\Delta}, \text{ with } W_{\Delta} = \frac{N}{T} \sum_{k=1}^{[T/\Delta]} (\bar{Z}_{k\Delta} - \bar{Z}_{(k-1)\Delta} - \Delta \hat{m})^2. $$

As it appears in the analysis below, $w_{\infty}^N$ is the limit of $W_{\Delta}$ when $T, \Delta \to \infty$. Hence, the choice $\hat{w} = W_{\Delta}$ would give a consistent estimator. However, the two estimators $W_{2\Delta}$ and $W_{\Delta}$ share a common bias which can be eliminated by considering $\hat{w} = 2W_{2\Delta} - W_{\Delta}$, which in turn drastically improves the rate of convergence. This common bias is related with the quantity $S_{\theta,t}^N$ defined in Lemma E.2.

Note that by the triangle inequality, the following inequality holds:

$$ (59) \qquad \left| \hat{w} - w_{\infty}^{N} \right| \le 2D_{2\Delta,T}^{N,1} + D_{\Delta,T}^{N,1} + 2D_{2\Delta,T}^{N,2} + D_{\Delta,T}^{N,2} + D_{T}^{N,3}, $$

where

$$ (60) \quad D_{\Delta,T}^{N,1} = \frac{N}{T} \left| \sum_{s=1}^{[T/\Delta]} (\bar{Z}_{s\Delta} - \bar{Z}_{(s-1)\Delta} - \Delta \hat{m})^2 - \sum_{s=1}^{[T/\Delta]} (\bar{Z}_{s\Delta} - \bar{Z}_{(s-1)\Delta} - \mathbb{E}_{\theta} [\bar{Z}_{s\Delta} - \bar{Z}_{(s-1)\Delta}])^2 \right|, $$

$$ (61) \quad D_{T}^{N,2} = \frac{N}{T} \left| \sum_{s=1}^{[T/\Delta]} (\bar{Z}_{s\Delta} - \bar{Z}_{(s-1)\Delta} - \mathbb{E}_{\theta} [\bar{Z}_{s\Delta} - \bar{Z}_{(s-1)\Delta}])^2 - \mathbb{E}_{\theta} \left[ \sum_{s=1}^{[T/\Delta]} (\bar{Z}_{s\Delta} - \bar{Z}_{(s-1)\Delta} - \mathbb{E}_{\theta} [\bar{Z}_{s\Delta} - \bar{Z}_{(s-1)\Delta}])^2 \right] \right|, $$

and

$$ (62) \quad D_{T}^{N,3} = \left| \frac{2N}{T} \mathbb{E}_{\theta} \left[ \sum_{s=1}^{[T/2\Delta]} (\bar{Z}_{2s\Delta} - \bar{Z}_{2(s-1)\Delta} - \mathbb{E}_{\theta} [\bar{Z}_{2s\Delta} - \bar{Z}_{2(s-1)\Delta}])^2 \right] - \frac{N}{T} \mathbb{E}_{\theta} \left[ \sum_{s=1}^{[T/\Delta]} (\bar{Z}_{s\Delta} - \bar{Z}_{(s-1)\Delta} - \mathbb{E}_{\theta} [\bar{Z}_{s\Delta} - \bar{Z}_{(s-1)\Delta}])^2 \right] - w_{\infty}^{N} \right|. $$

In the rest of this section, propositions provide upper-bounds for the terms involved in the right-hand side of (59), while lemmas provide intermediate results.

**PROPOSITION E.1.** There exists a constant $K$ depending only on $\lambda$ such that for all $N \ge 1$, $T \ge 2$ and $1 \le \Delta \le [T/2]$, we have that

$$ \mathbb{E}_{\theta} [D_{2\Delta,T}^{N,1} + D_{\Delta,T}^{N,1}] \le K \frac{\Delta}{T}. $$
PROOF. Using that $\sum_{s=1}^{|T/\Delta|} (\bar{Z}_{s\Delta} - \bar{Z}_{(s-1)\Delta}) = |T/\Delta|\Delta\hat{m}$ and $\mathbb{E}_{\theta}[\bar{Z}_{\Delta}] = \Delta\mathbb{E}_{\theta}[\hat{m}]$, one can check that

$$D_{\Delta,T}^{N,1} = \frac{N}{T} \Delta^2 \left| \frac{T}{\Delta} \right| |\hat{m} - \mathbb{E}_{\theta}[\hat{m}]|^2 \le N \Delta |\hat{m} - \mathbb{E}_{\theta}[\hat{m}]|^2,$$

which implies that $\mathbb{E}_{\theta}[D_{\Delta,T}^{N,1}] \le N\Delta \mathrm{Var}_{\theta}(\hat{m})$. Replacing $\Delta$ by $2\Delta$ in the previous identity, we obtain $\mathbb{E}_{\theta}[D_{2\Delta,T}^{N,1}] \le 2N\Delta \mathrm{Var}_{\theta}(\hat{m})$. Hence, the result follows from Equation (32). $\square$

To deal with the other two terms, we need to obtain a fine estimate on $\mathrm{Var}_{\theta}(\bar{Z}_{\Delta}) = \mathbb{E}_{\theta}[(\bar{U}_{\Delta})^2]$ (recall (23)). This can be done as follows. Recall that for all $t \ge 1$ and $1 \le i \le N$,

$$M_{i,t} = \sum_{s=1}^{t} (X_{i,s} - p_{\theta,i}(X_{s-1})), \text{ and } M_t = (M_{1,t}, \dots, M_{N,t}).$$

Let us denote $\mathcal{F}_t = \sigma(X_0, \dots, X_t)$. From the fact that the process $(X_t)_{t \in \mathbb{Z}}$ is Markovian, it follows that $(M_t)_{t \ge 1}$ is a martingale with respect to the filtration $(\mathcal{F}_t)_{t \ge 1}$.

Recall that $Q^N = (I_N - (1-\lambda)A^N)^{-1}$. Then (24) implies

$$(63) \qquad U_t = Q^N M_t + (Q^N - I_N)(X_0 - X_t), \quad t \ge 1.$$

LEMMA E.2. There exists a constant $K$ depending only on $\lambda$ such that for all $t \ge 1$ and $N \ge 1$,

$$\mathbb{E}_{\theta} \left[ \left( \overline{U}_t \right)^2 \right] = \mathbb{E}_{\theta} \left[ \left( \overline{\left( Q^N M_t \right)} \right)^2 \right] + S_{\theta,t}^N + r_{\theta,t}^N,$$

where $|r_{\theta,t}^N| \le K(1-\lambda)^t N^{-1}$ and $S_{\theta,t}^N$ has the following properties: $|S_{\theta,t}^N| \le KN^{-1}$ and $|S_{\theta,2t}^N - S_{\theta,t}^N| \le K(1-\lambda)^t N^{-1}$.

PROOF. Starting from Equation (63), one can check that for any $t \ge 1$,

$$\mathbb{E}_{\theta} \left[ \left( \overline{U}_t \right)^2 \right] = \mathbb{E}_{\theta} \left[ \left( \overline{\left( Q^N M_t \right)} \right)^2 \right] + H_{\theta,t}^N,$$

where

$$H_{\theta,t}^N = 2 \mathbb{E}_{\theta} \left[ \overline{Q^N M_t} \overline{(Q^N - I_N)(X_0 - X_t)} \right] + \mathbb{E}_{\theta} \left[ \left( \overline{(Q^N - I_N)(X_0 - X_t)} \right)^2 \right].$$

Hence, to conclude the proof, it remains to show that we can write $H_{\theta,t}^N = S_{\theta,t}^N + r_{\theta,t}^N$ where $|r_{\theta,t}^N| \le K(1-\lambda)^t N^{-1}$ and $S_{\theta,t}^N$ is such that $|S_{\theta,2t}^N - S_{\theta,t}^N| \le K(1-\lambda)^t N^{-1}$, for all $t \ge 1$ and $N \ge 1$. To see that, first of all, observe that

$$\begin{aligned} \overline{(Q^N - I_N)(X_0 - X_t)} &= N^{-1} \sum_{i=1}^{N} \sum_{j=1}^{N} (Q^N(i,j) - I_N(i,j))(X_{j,0} - X_{j,t}) \\ &= N^{-1} \sum_{j=1}^{N} (X_{j,0} - X_{j,t})(c_j^N - 1), \end{aligned}$$
where $c^N = (Q^N)^{\top}1_N$, so that

$$ \begin{aligned} \overline{(Q^N - I_N)(X_0 - X_t)} \bigg|^2 &= N^{-2} \sum_{i=1}^{N} \sum_{j=1}^{N} (c_i^N - 1)(c_j^N - 1)(X_{i,0} - X_{i,t})(X_{j,0} - X_{j,t}) \\ &= N^{-2} \sum_{i=1}^{N} \sum_{j=1}^{N} (c_i^N - 1)(c_j^N - 1)(Y_{i,0} - Y_{i,t})(Y_{j,0} - Y_{j,t}), \end{aligned} $$

where the second equality follows from the definition of $Y_{i,t} = X_{i,t} - \mathbb{E}_{\theta}[X_{i,t}]$ and the stationarity of the process $(X_t)_{t \in \mathbb{Z}}$. Using the stationarity once more and the above equation, we then deduce that

$$ \mathbb{E}_{\theta} \left[ \left( \overline{(Q^N - I_N)(X_0 - X_t)} \right)^2 \right] = 2N^{-2} \sum_{i=1}^{N} \sum_{j=1}^{N} (c_i^N - 1)(c_j^N - 1) \left[ \begin{aligned} &\mathrm{Cov}_{\theta}(Y_{i,0}, Y_{j,0}) \\ &- \mathrm{Cov}_{\theta}(Y_{i,0}, Y_{j,t}) \end{aligned} \right], $$

so that by applying Lemma 3.7 and using Inequality (84), we obtain that

$$ \mathbb{E}_{\theta} \left[ \left( \overline{(Q^N - I_N)(X_0 - X_t)} \right)^2 \right] = \frac{2}{N^2} \sum_{i=1}^{N} \sum_{j=1}^{N} (c_i^N - 1)(c_j^N - 1) \mathrm{Cov}_{\theta}(Y_{i,0}, Y_{j,0}) + r_{\theta,t}^N, $$

where $|r_{\theta,t}^N| \le K(1 - \lambda)^t N^{-1}$.

Now, proceeding similarly as above, one can check that

$$ \mathbb{E}_{\theta} \left[ \overline{Q^N M_t (Q^N - I_N)(X_0 - X_t)} \right] = N^{-2} \sum_{i=1}^{N} \sum_{j=1}^{N} c_i^N (c_j^N - 1) \mathbb{E}_{\theta} [M_{i,t}(X_{j,0} - X_{j,t})]. $$

For each $1 \le i, j \le N$ and $t \ge 1$, denote $S_{\theta,t}^N(i, j) = \mathbb{E}_{\theta}[M_{i,t}(X_{j,0} - X_{j,t})]$. We claim that there exists a constant $K > 0$ such that for all $1 \le i, j \le N$ and $t \ge 1$, the following inequality holds: $|S_{\theta,t}^N(i, j)| \le K(\delta_{ij} + N^{-1})$ and

$$ |S_{\theta,2t}^N(i,j) - S_{\theta,t}^N(i,j)| \le K(1-\lambda)^t N^{-1}. $$

Once this claim is proved, using Inequality (84) and Lemma 3.7 it is immediate to check that $H_{\theta,t}^N = S_{\theta,t}^N + r_{\theta,t}^N$, where

$$ S_{\theta,t}^N = N^{-2} \sum_{i=1}^{N} \sum_{j=1}^{N} [c_i^N(c_j^N - 1)S_{\theta,t}^N(i,j) + 2(c_i^N - 1)(c_j^N - 1)\mathrm{Cov}_{\theta}(X_{i,0}, X_{j,0})] $$

satisfies $|S_{\theta,2t}^N - S_{\theta,t}^N| \le K(1 - \lambda)^t N^{-1}$ and $|S_{\theta,t}^N| \le KN^{-1}$, implying the result.

To prove the claim above, we first write

$$ S_{\theta,t}^N(i,j) = \mathbb{E}_{\theta} [M_{i,t}(X_{j,0} - X_{j,t})] = \sum_{s=1}^{t} \mathbb{E}_{\theta} [(X_{i,s} - p_{\theta,i}(X_{s-1}))(X_{j,0} - X_{j,t})] $$

and then use the fact that

$$ \mathbb{E}_{\theta}[(X_{i,s} - p_{\theta,i}(X_{s-1}))X_{j,0}] = \mathbb{E}_{\theta}[X_{j,0}\mathbb{E}_{\theta}[(X_{i,s} - p_{\theta,i}(X_{s-1}))|\mathcal{F}_{s-1}]] = 0 $$

to deduce that

$$ \mathbb{E}_{\theta} [M_{i,t}(X_{j,0} - X_{j,t})] = \sum_{s=1}^{t} (\mathbb{E}_{\theta} [p_{\theta,i}(X_{s-1})X_{j,t}] - \mathbb{E}_{\theta} [X_{i,s}X_{j,t}]). $$
Next, we use (1) and the stationarity to obtain that

$$ \sum_{s=1}^{t} \mathbb{E}_{\theta} [p_{\theta,i}(X_{s-1})X_{j,t}] = tm_j^N(\mu - (1-\lambda)L_i^{N,\bullet-}) \\ + (1-\lambda) \sum_{k=1}^{N} A^N(i,k) \sum_{s=1}^{t} \mathbb{E}_{\theta} [X_{k,s-1}X_{j,t}]. $$

Hence, using that $\mathbb{E}_{\theta}[X_{k,s-1}X_{j,t}] = \text{Cov}_{\theta}(X_{k,s-1}, X_{j,t}) + m_k^N m_j^N$ and the stationarity once more, one can check that

$$ \sum_{s=1}^{t} \mathbb{E}_{\theta} [p_{\theta,i}(X_{s-1})X_{j,t}] = tm_j^N(\mu - (1-\lambda)L_i^{N,\bullet-} + (1-\lambda)(A^N m^N)_i) \\ + (1-\lambda) \sum_{k=1}^{N} A^N(i,k) \sum_{s=1}^{t} \text{Cov}_{\theta}(X_{k,0}, X_{j,s}), $$

and similarly that

$$ \sum_{s=1}^{t} \mathbb{E}_{\theta} [X_{i,s}X_{j,t}] = tm_i^N m_j^N - \sum_{s=0}^{t-1} \text{Cov}_{\theta} [X_{i,0}, X_{j,s}]. $$

Therefore, combining the last two equations and using (14), we obtain that

$$ S_{\theta,t}^{N}(i,j) = (1-\lambda) \sum_{k=1}^{N} A^{N}(i,k) \sum_{s=1}^{t} \text{Cov}_{\theta}(X_{k,0}, X_{j,s}) + \sum_{s=0}^{t-1} \text{Cov}_{\theta}(X_{i,0}, X_{j,s}), $$

so that Lemma 3.7 allows us to conclude that

$$ \begin{aligned} |S_{\theta,2t}^N(i,j) - S_{\theta,t}^N(i,j)| & \le \\ & \le (1-\lambda) \sum_{k=1}^N |A^N(i,k)| \sum_{s=t+1}^{2t} |\text{Cov}_{\theta}(X_{k,0}, X_{j,s})| + \sum_{s=t}^{2t-1} |\text{Cov}_{\theta}(X_{i,0}, X_{j,s})| \\ & \le \frac{K(1-\lambda)^t}{N} \left( \sum_{s=0}^{t-1} (1-\lambda)^s \right) [(1-\lambda)^2 N \|A^N\|_{\infty} + 1] \le \frac{K(1-\lambda)^t}{N}. \end{aligned} $$

Finally, using once more Lemma 3.7, one can easily check that $|S_{\theta,t}^N(i,j)| \le K(\delta_{ij} + N^{-1})$ proving the claim. □

The next result will be important to find the leading term of $\mathbb{E}_{\theta} \left[ \left( \overline{Q^N M_t} \right)^2 \right]$.

LEMMA E.3. For all $t \ge 1$, the following equality holds.

$$ \mathbb{E}_{\theta} \left[ (X_t - p_{\theta}(X_{t-1}))^2 \right] = m^N - (m^N)^2 + r_{\theta}^N, $$
where $\|r_{\theta}^N\|_{\infty} \le KN^{-1}$ for some constant $K > 0$ depending only on $\lambda$.

REMARK E.4. Note that $X_t$, $p_{\theta}(X_{t-1})$ and $m^N$ are vectors so that the square terms in Lemma E.3 have to be interpreted in terms of the Hadamard product of vectors. Also, the remainder term $r_{\theta}^N$ does not depend on time $t$ because the process $(X_t)_{t \in \mathbb{Z}}$ is assumed to be stationary.
PROOF. First, observe that

$$
\begin{align*}
\mathbb{E}_{\theta} \left[ (X_{i,t} - p_{i,\theta}(X_{t-1}))^2 \right] &= \mathbb{E}_{\theta} [X_{i,t}] - \mathbb{E}_{\theta} [p_{i,\theta}^2(X_{t-1})] \\
&= m_i^N - \mathbb{E}_{\theta} [p_{i,\theta}^2(X_0)],
\end{align*}
$$

where the second equality follows by the stationarity. Now, from (1), we have that

$$
p_{i,\theta}(X_0) = (\mu - (1 - \lambda)L_i^{N, \bullet -}) + (1 - \lambda)(A^N X_0)_i,
$$

so that

$$
(64) \quad \begin{aligned}
\mathbb{E}_{\theta} [p_{i,\theta}^2(X_0)] = (\mu - (1-\lambda)L_i^{N,\bullet-})^2 &+ 2(\mu - (1-\lambda)L_i^{N,\bullet-})(1-\lambda)(A^N \mathbb{E}_{\theta}[X_0])_i \\
&+ (1-\lambda)^2 \mathbb{E}_{\theta} [((A^N X_0)_i)^2].
\end{aligned}
$$

Next, notice that

$$
\mathbb{E}_{\theta} [((A_{\theta}^{N} X_{0})_{i})^{2}] = \sum_{j=1}^{N} \sum_{k=1}^{N} A^{N}(i,j)A^{N}(i,k) \mathbb{E}_{\theta} [X_{k,0}X_{j,0}],
$$

so that using that $\mathbb{E}_{\theta}[X_{k,0}X_{j,0}] = m_j^N m_k^N + \text{Cov}_{\theta}(X_{k,0}, X_{j,0})$, we deduce that

$$
(65) \quad \mathbb{E}_{\theta} [((A^N X_0)_i)^2] = \left( \sum_{j=1}^{N} A^N(i,j) m_j^N \right)^2 + \tilde{r}_{\theta,i}^N
$$

where

$$
\tilde{r}_{\theta,i}^{N} = \sum_{j=1}^{N} \sum_{k=1}^{N} A^{N}(i,j)A^{N}(i,k) \mathrm{Cov}_{\theta}(X_{k,0}, X_{j,0}).
$$

Using the fact that $|A^N(i,j)| \le 1/N$ and Lemma 3.7, we have

$$
|\tilde{r}_{\theta,i}^N| \le N^{-2} \sum_{j=1}^N \sum_{k=1}^N |\mathrm{Cov}_\theta[X_{k,0}, X_{j,0}]| \le N^{-2}(N + N^2 K(1-\lambda)N^{-1}) \le KN^{-1}.
$$

Finally, recalling that $m^N = \mathbb{E}_{\theta}[X_0]$, it follows from (64) and (65) that

$$
\mathbb{E}_{\theta} [p_{i,\theta}^2(X_0)] = \left[ (\mu - (1-\lambda)L_i^{N,\bullet-}) + (1-\lambda)(A^N m^N)_i \right]^2 + (1-\lambda)^2 \tilde{r}_{\theta,i}^N.
$$

The result then follows from equation (15).

<figure><img src="image_13.png" /></figure>

As a consequence of Lemma E.3, we will now find the leading term of $\mathbb{E}_{\theta} \left[ \left( \overline{Q^N M_t} \right)^2 \right]$.

LEMMA E.5. For any $t \ge 1$ and $N \ge 1$,

$$
\mathbb{E}_{\theta} \left[ \left| \overline{Q^N M_t} \right|^2 \right] = \frac{t}{N} w_{\infty}^N + r_{\theta,t}^N,
$$

where $w_{\infty}^N$ is defined in (31) and $r_{\theta,t}^N$ is such that $|r_{\theta,t}^N| \le KtN^{-2}$ for some $K$ depending only on $\lambda$.
PROOF. We start observing that

$$ \overline{Q^N M_t} = N^{-1} \sum_{i=1}^{N} \sum_{j=1}^{N} Q^N(i,j) M_{j,t} = N^{-1} \sum_{j=1}^{N} c_j^N M_{j,t}, $$

and then we use that for all $s \neq h$ and $1 \le j, k \le N$,

$$ \mathbb{E}_{\theta}[(X_{i,s} - p_{i,\theta}(X_{s-1}))(X_{k,h} - p_{k,\theta}(X_{h-1}))] = 0, $$

to deduce that

$$ \begin{aligned} \mathbb{E}_{\theta} \left[ \left| \overline{Q^N M_t} \right|^2 \right] &= N^{-2} \sum_{j=1}^{N} \sum_{k=1}^{N} c_j^N c_k^N \mathbb{E}_{\theta} [M_{j,t} M_{k,t}] \\ &= N^{-2} \sum_{j=1}^{N} (c_j^N)^2 \sum_{s=1}^{t} \mathbb{E}_{\theta} \left[ (X_{j,s} - p_{j,\theta}(X_{s-1}))^2 \right]. \end{aligned} $$

Hence, it follows from Lemma E.3 that

$$ \begin{aligned} \mathbb{E}_{\theta} \left[ \left| \overline{Q^N M_t} \right|^2 \right] &= N^{-2} \sum_{j=1}^{N} (c_j^N)^2 \sum_{s=1}^{t} \mathbb{E}_{\theta} \left[ (X_{j,s} - p_{j,\theta}(X_{s-1}))^2 \right] \\ &= tN^{-2} \sum_{j=1}^{N} (c_j^N)^2 m_j^N (1 - m_j^N) + tN^{-2} \langle (c^N)^2, r_{\theta}^N \rangle \\ &= tN^{-1} w_{\infty}^N + tN^{-2} \langle (c^N)^2, r_{\theta}^N \rangle. \end{aligned} $$

Finally, by applying Hölder inequality and then using that $\|v\|_1 \le N\|v\|_{\infty}$ for $v \in \mathbb{R}^N$, we obtain

$$ \langle (c^N)^2, r_{\theta}^N \rangle \le \|(c^N)^2\|_{\infty} \|r_{\theta}^N\|_1 \le \|(c^N)^2\|_{\infty} N \|r_{\theta}^N\|_{\infty}, $$

and the result follows from inequality (84) and Lemma E.3. □

For later use, let us mention the following immediate corollary of Lemmas E.2 and E.5.

LEMMA E.6. There exists a constant $K > 0$ depending only on $\lambda$ such that for all $t \ge 1$ and $N \ge 1$,

$$ \mathbb{E}_{\theta} \left[ \left( \overline{U}_t \right)^2 \right] \le K t N^{-1}. $$

PROOF. First, observe that $w_{\infty}^N$ defined in (31) satisfies $|w_{\infty}^N| \le \lambda^{-2}/4$ by Inequality (84), and then combine this fact with Lemmas E.2 and E.5. □

We are now in position to deal with the term $D_T^{N,3}$ defined in (62).

PROPOSITION E.7. There exists a constant $K > 0$ depending only on $\lambda$ such that for all $N \ge 1$, $T \ge 2$ and $1 \le \Delta \le \lfloor T/2 \rfloor$,

$$ \mathbb{E}_{\theta} \left[ D_T^{N,3} \right] \le K \left( \frac{1}{N} + \frac{(1-\lambda)^{\Delta}}{\Delta} + \frac{\Delta}{T} \right). $$
PROOF. First, observe that by stationarity,

$$D_T^{N,3} = \left| \frac{2N}{T} \lfloor T/2\Delta \rfloor \mathbb{E}_{\theta} \left[ (\bar{U}_{2\Delta})^2 \right] - \frac{N}{T} \lfloor T/\Delta \rfloor \mathbb{E}_{\theta} \left[ (\bar{U}_{\Delta})^2 \right] - w_{\infty}^N \right|.$$

Then, by using the triangle inequality, one can check that

$$\left| \frac{2N}{T} \lfloor T/2\Delta \rfloor \mathbb{E}_{\theta} \left[ (\bar{U}_{2\Delta})^2 \right] - \frac{N}{T} \lfloor T/\Delta \rfloor \mathbb{E}_{\theta} \left[ (\bar{U}_{\Delta})^2 \right] - w_{\infty}^N \right| \\ \leq \left| \frac{N}{\Delta} \left( \mathbb{E}_{\theta} \left[ (\bar{U}_{2\Delta})^2 \right] - \mathbb{E}_{\theta} \left[ (\bar{U}_{\Delta})^2 \right] \right) - w_{\infty}^N \right| + \frac{2N}{T} \left( \mathbb{E}_{\theta} \left[ (\bar{U}_{2\Delta})^2 \right] + \mathbb{E}_{\theta} \left[ (\bar{U}_{\Delta})^2 \right] \right).$$

Then Lemmas E.2, E.5 and E.6 imply the result. □

It remains to deal with $D_{\Delta,T}^{N,2}$ and $D_{2\Delta,T}^{N,2}$. To that end, we shall use two additional lemmas. The first one is the following.

LEMMA E.8. There exists a constant $K > 0$ depending only on $\lambda$ such that for all $N \ge 1$, $\Delta \ge 1$ and $t \ge 2$,

$$\mathrm{Cov}_{\theta} \left[ (\bar{U}_{\Delta})^2, (\bar{U}_{t\Delta} - \bar{U}_{(t-1)\Delta})^2 \right] \leq K N^{-2} (1 - \lambda)^{(t-2)\Delta}.$$

PROOF. First, use that $\bar{U}_t = N^{-1} \sum_{i=1}^N \sum_{s=1}^t Y_{i,s}$ to deduce that

$$\mathrm{Cov}_{\theta} \left[ (\bar{U}_{\Delta})^2, (\bar{U}_{t\Delta} - \bar{U}_{(t-1)\Delta})^2 \right] = N^{-4} \sum_{z_1 \in F_1} \sum_{z_2 \in F_1} \sum_{z_3 \in F_t} \sum_{z_4 \in F_t} \mathrm{Cov}_{\theta} [Y_{z_1} Y_{z_2}, Y_{z_3} Y_{z_4}],$$

where $z_k = (i_k, t_k)$ for $k \in \{1, \dots, 4\}$, and $F_s := \{1, \dots, N\} \times \{(s-1)\Delta + 1, \dots, s\Delta\}$ for $s \ge 1$. In the rest of the proof, we shall denote $\mathbf{z} = (z_1, z_2, z_3, z_4)$ and $B(\mathbf{z}) = \mathrm{Cov}_{\theta} [Y_{z_1} Y_{z_2}, Y_{z_3} Y_{z_4}]$. Notice that the covariance function $B(\mathbf{z})$ is symmetric with respect to $z_1$ and $z_2$, as well as with respect to $z_3$ and $z_4$. For each $p \in \{1, \dots, 6\}$ and $s \ge 1$, we denote $C_{p,s}$ the set of all vectors $\mathbf{z} \in F_1 \times F_1 \times F_s \times F_s$ satisfying the conditions of Item $p$ of Lemma 3.8 and such that $t_1 \le t_2$ and $t_3 \le t_4$. With this notation, one can check using the symmetry of the covariance function that the following inequality holds:

$$(66) \qquad \mathrm{Cov}_{\theta} \left[ (\bar{U}_{\Delta})^2, (\bar{U}_{t\Delta} - \bar{U}_{(t-1)\Delta})^2 \right] \le 4N^{-4} \sum_{p=1}^{6} \sum_{\mathbf{z} \in C_{p,t}} |B(\mathbf{z})|.$$

Notice that $C_{1,t} = C_{2,t} = C_{5,t} = \emptyset$ for $t \ge 2$, so that we need only to control the sum of covariances for the cases $p \in \{3, 4, 6\}$. This is done in the 3 steps below.

Step 1. Here, we show that there exists a constant $K > 0$ such that for all $N \ge 1$, $\Delta \ge 1$ and $t \ge 2$,

$$\sum_{\mathbf{z} \in C_{3,t}} |B(\mathbf{z})| \le KN(1-\lambda)^{(t-2)\Delta+1}.$$

First, note that Item 3 of Lemma 3.8 implies that

$$\sum_{\mathbf{z} \in C_{3,t}} |B(\mathbf{z})| \le KN^2 \sum_{t_1=1}^{\Delta} \sum_{t_3=(t-1)\Delta+1}^{t\Delta} (1-\lambda)^{(t_3-t_1)} N^{-1} \\ = KN(1-\lambda)^{(t-2)\Delta+1} \sum_{k=0}^{\Delta-1} (1-\lambda)^k \sum_{t_1=1}^{\Delta} (1-\lambda)^{\Delta-t_1}.$$
Step 1 follows by noticing that $\sum_{t_1=1}^{\Delta} (1-\lambda)^{\Delta-t_1} = \sum_{k=0}^{\Delta-1} (1-\lambda)^k$, and then by using that $\sum_{k=0}^{\Delta-1} (1-\lambda)^k \le \sum_{k=0}^{\infty} (1-\lambda)^k \le \lambda^{-1}$.

Step 2. Here, we prove that there exists a constant $K > 0$ such that for all $N \ge 1$, $\Delta \ge 1$ and $t \ge 2$,

$$ \sum_{\mathbf{z} \in C_{4,t}} |B(\mathbf{z})| \le KN(1-\lambda)^{(t-2)\Delta+1}. $$

First, by Remark 3.9, we can write

$$ \sum_{\mathbf{z} \in C_{4,t}} |B(\mathbf{z})| = \sum_{\mathbf{z} \in C_{4,t}} \mathbf{1}_{\{z_1=z_2, z_3 \neq z_4\}} |B(\mathbf{z})| + \sum_{\mathbf{z} \in C_{4,t}} \mathbf{1}_{\{z_1 \neq z_2, z_3=z_4\}} |B(\mathbf{z})|. $$

To conclude, we apply Item 4 of Lemma 3.8 to obtain an upper bound for each one of the terms on the right hand side of the above identity. Since both terms are treated very similarly, we explain how we handle the first term only. By applying Item 4 of Lemma 3.8 to the first one, we obtain that

$$ \sum_{\mathbf{z} \in C_{4,t}} \mathbf{1}_{\{z_1=z_2, z_3 \neq z_4\}} |B(\mathbf{z})| \le KN^3 \sum_{t_1=1}^{\Delta} \sum_{t_3=(t-1)\Delta+1}^{t\Delta} \sum_{t_4=t_3}^{t\Delta} (1-\lambda)^{(t_4-t_3)+(t_3-t_1)N^{-2}}. $$

Now, since for any $(t-1)\Delta + 1 \le t_3 \le t\Delta$,

$$ (67) \quad \sum_{t_4=t_3}^{t\Delta} (1-\lambda)^{(t_4-t_3)} = \sum_{k=0}^{t\Delta-t_3} (1-\lambda)^k \le \sum_{k=0}^{\Delta-1} (1-\lambda)^k \le \lambda^{-1}, $$

we can then proceed as in Step 1 to conclude that

$$ \sum_{\mathbf{z} \in C_{4,t}} \mathbf{1}_{\{z_1=z_2, z_3 \neq z_4\}} |B(\mathbf{z})| \le KN(1-\lambda)^{(t-2)\Delta+1}. $$

Step 3. Here, we prove that there exists a constant $K > 0$ such that for all $N \ge 1$, $\Delta \ge 1$ and $t \ge 2$,

$$ \sum_{\mathbf{z} \in C_{6,t}} |B(\mathbf{z})| \le KN^2(1-\lambda)^{(t-2)\Delta+1}. $$

By applying Item 6 of Lemma 3.8 and inequality (67), one can check that

$$ \sum_{\mathbf{z} \in C_{6,t}} |B(\mathbf{z})| \le KN^4 \sum_{t_1=1}^{\Delta} \sum_{t_2=t_1}^{\Delta} \sum_{t_3=(t-1)\Delta+1}^{t\Delta} N^{-2}(1-\lambda)^{t_3-t_1}. $$

Some algebraic computations imply that

$$ \sum_{t_1=1}^{\Delta} \sum_{t_2=t_1}^{\Delta} \sum_{t_3=(t-1)\Delta+1}^{t\Delta} (1-\lambda)^{t_3-t_1} \le (1-\lambda)^{(t-2)\Delta+1} \left( \sum_{k=0}^{\Delta-1} (1-\lambda)^k \right) \times \left( \sum_{k=0}^{\Delta-1} (k+1)(1-\lambda)^k \right). $$

Putting together these estimates and using the fact that $\sum_{k=0}^{\infty} (k+1)(1-\lambda)^k < \infty$, we conclude the proof of Step 3.

Finally, combining Steps 1 through 3 and (66), the result follows. $\square$

<figure><img src="image_14.png" /></figure>
LEMMA E.9. There exists a constant K depending only on $\lambda, \mu$ and $p$, such that for all $N \ge 1$ and $\Delta \ge 1$,

$$ \mathrm{Var}_{\theta}[(\bar{U}_{\Delta})^2] \le K\Delta^2 N^{-2}. $$

PROOF. We use the notation of the proof of Lemma E.8, except that we now denote $C_{1,6}$ the set of all vectors $\mathbf{z} \in \mathbf{F}_1 \times \mathbf{F}_1 \times \mathbf{F}_1 \times \mathbf{F}_1$ satisfying the conditions of Item 7. of Lemma 3.8, that is, without imposing the constraint that $t_1 > t_4$ or $t_3 > t_2$. By taking $t = 1$ in inequality (66), it follows that

$$ (68) \qquad \mathrm{Var}_{\theta} \left[ (\bar{U}_{\Delta})^2 \right] \le 4N^{-4} \left( \sum_{p=1}^{6} \sum_{\mathbf{z} \in C_{1,p}} |B(\mathbf{z})| \right). $$

So we need to show that the sum of covariances is at most $K\Delta^2 N^2$. This is done in the 5 steps below.

Step 1. First, we show that there exists a constant $K > 0$ such that for all $N \ge 1$ and $\Delta \ge 1$,

$$ \sum_{p=1}^{2} \sum_{\mathbf{z} \in C_{1,p}} |B(\mathbf{z})| \le K(\Delta N)^2. $$

By Item 1 and Item 2 of Lemma 3.8, we have that

$$ \sum_{p=1}^{2} \sum_{\mathbf{z} \in C_{1,p}} |B(\mathbf{z})| \le |C_{1,1}| + |C_{1,2}| \le \Delta N + 2\Delta N(\Delta N - 1) \le 3(\Delta N)^2, $$

so that the result follows with $K = 3$.

Step 2. Next, we show that there exists a constant $K > 0$ such that for all $N \ge 1$ and $\Delta \ge 1$,

$$ \sum_{\mathbf{z} \in C_{1,3}} |B(\mathbf{z})| \le KN\Delta^2. $$

Item 3 of Lemma 3.8, ensures that

$$ \sum_{\mathbf{z} \in C_{1,3}} |B(\mathbf{z})| \le KN^2 \sum_{t_1=1}^{\Delta} \sum_{t_3=1}^{\Delta} (1-\lambda)^{|t_3-t_1|} N^{-1}. $$

Since $\sum_{t_1=1}^{\Delta} \sum_{t_3=1}^{\Delta} (1-\lambda)^{|t_3-t_1|} \le \Delta^2$, the result follows from the previous inequality.

Step 3. Here, we show that there exists a constant $K > 0$ such that for all $N \ge 1$ and $\Delta \ge 1$,

$$ \sum_{\mathbf{z} \in C_{1,4}} |B(\mathbf{z})| \le KN\Delta^2. $$

We proceed very similarly as in the proof of Step 2 of Lemma E.8. First, observe that we can write (see Remark 3.9),

$$ \sum_{\mathbf{z} \in C_{1,4}} |B(\mathbf{z})| = \sum_{\mathbf{z} \in C_{1,4}} \mathbf{1}_{\{z_1=z_2, z_3 \neq z_4\}} |B(\mathbf{z})| + \sum_{\mathbf{z} \in C_{1,4}} \mathbf{1}_{\{z_1 \neq z_2, z_3=z_4\}} |B(\mathbf{z})|. $$
To obtain the result, we apply Item 4 of Lemma 3.8 to obtain an upper bound for each term on the right hand side of the above identity. Since both terms are treated very similarly, we explain how we deal with the first term only. By applying Point 4 of Lemma 3.8 to the first one, we obtain that

$$ \sum_{\mathbf{z} \in C_{1,4}} \mathbf{1}_{\{z_1=z_2, z_3 \neq z_4\}} |B(\mathbf{z})| \leq K N^3 \sum_{t_1=1}^{\Delta} \sum_{t_3=1}^{\Delta} \sum_{t_4=1}^{\Delta} \mathbf{1}_{\{t_3 \leq t_4\}} (1-\lambda)^{|s_3-s_2|+|s_2-s_1|} N^{-2}, $$

where $s_1 \leq s_2 \leq s_3$ denotes the ordering of the triple $(t_1, t_3, t_4)$. Next, we use that $|s_3 - s_2| + |s_2 - s_1| \geq t_4 - t_3$ to deduce that

$$ \begin{aligned} \sum_{t_1=1}^{\Delta} \sum_{t_3=1}^{\Delta} \sum_{t_4=1}^{\Delta} \mathbf{1}_{\{t_3 \le t_4\}} (1-\lambda)^{|s_3-s_2|+|s_2-s_1|} &\le \Delta \sum_{t_3=1}^{\Delta} \sum_{t_4=1}^{\Delta} \mathbf{1}_{\{t_3 \le t_4\}} (1-\lambda)^{t_4-t_3} \\ &\le \Delta^2 \sum_{k=0}^{\Delta-1} (1-\lambda)^k \le \lambda^{-1} \Delta^2. \end{aligned} $$

As a consequence, we obtain that

$$ \sum_{\mathbf{z} \in C_{1,4}} \mathbf{1}_{\{z_1=z_2, z_3 \neq z_4\}} |B(\mathbf{z})| \leq K N \lambda^{-1} \Delta^2, $$

and the result follows.

Step 4. In this step, we show that there exists a constant $K > 0$ such that for all $N \geq 1$ and $\Delta \geq 1$,

$$ \sum_{\mathbf{z} \in C_{1,5}} |B(\mathbf{z})| \leq K N^2 \Delta^2. $$

We start observing that

$$ \sum_{\mathbf{z} \in C_{1,5}} |B(\mathbf{z})| = \sum_{\mathbf{z} \in C_{1,5}} \mathbf{1}_{\{z_1=z_3, z_2 \neq z_4\}} |B(\mathbf{z})| + \sum_{\mathbf{z} \in C_{1,5}} \mathbf{1}_{\{z_1 \neq z_3, z_2 = z_4\}} |B(\mathbf{z})|. $$

Hence, it suffices to provide an upper bound for each term on the right hand side of the above inequality. We will explain how we deal with the first one only. The second one can be treated similarly. By Item 5 of Lemma 3.8, we have that

$$ \sum_{\mathbf{z} \in C_{1,5}} \mathbf{1}_{\{z_1=z_3, z_2 \neq z_4\}} |B(\mathbf{z})| \leq K N^3 \sum_{t_1=1}^{\Delta} \sum_{t_2=1}^{\Delta} \sum_{t_4=1}^{\Delta} \mathbf{1}_{\{t_2 \leq t_4\}} (1-\lambda)^{(t_4-t_2)} N^{-1}. $$

Hence, by observing that

$$ \sum_{t_1=1}^{\Delta} \sum_{t_2=1}^{\Delta} \sum_{t_4=1}^{\Delta} \mathbf{1}_{\{t_2 \le t_4\}} (1-\lambda)^{(t_4-t_2)} \le \Delta^2 \sum_{k=0}^{\Delta-1} (1-\lambda)^k \le \Delta^2 \lambda^{-1}, $$

the result follows.

Step 5. Finally, we show that there exists a constant $K > 0$ such that for all $N \geq 1$ and $\Delta \geq 1$,

$$ \sum_{\mathbf{z} \in C_{1,6}} |B(\mathbf{z})| \leq K N^2 \Delta^2. $$
We start from

$$ \sum_{\mathbf{z} \in C_{1,6}} |B(\mathbf{z})| \le \sum_{\mathbf{z} \in C_{1,6}} \left[ \mathbf{1}_{\{t_3 \le t_4 \le t_1 \le t_2\}} + \mathbf{1}_{\{t_1 \le t_2 \le t_3 \le t_4\}} \right] |B(\mathbf{z})| + \sum_{\mathbf{z} \in C_{1,6}} \left[ \mathbf{1}_{\{t_1 \le t_3 \le t_2 \le t_4\}} + \mathbf{1}_{\{t_3 \le t_1 \le t_4 \le t_2\}} + \mathbf{1}_{\{t_1 \le t_3 \le t_4 \le t_2\}} + \mathbf{1}_{\{t_3 \le t_1 \le t_2 \le t_4\}} \right] |B(\mathbf{z})|. $$

In the sequel, we handle the two terms on the right-hand side of the above inequality separately.

Step 5.1 Here, we show that there exists a constant $K > 0$ such that for all $N \ge 1$ and $\Delta \ge 1$,

$$ \sum_{\mathbf{z} \in C_{1,6}} \left[ \mathbf{1}_{\{t_3 \le t_4 \le t_1 \le t_2\}} + \mathbf{1}_{\{t_1 \le t_2 \le t_3 \le t_4\}} \right] |B(\mathbf{z})| \le K N^2 \Delta. $$

By Item 6 of Lemma 3.8,

$$ \sum_{\mathbf{z} \in C_{1,6}} \mathbf{1}_{\{t_3 \le t_4 \le t_1 \le t_2\}} |B(\mathbf{z})| \le K N^4 \sum_{t_1=1}^{\Delta} \sum_{t_2=1}^{\Delta} \sum_{t_3=1}^{\Delta} \sum_{t_4=1}^{\Delta} \mathbf{1}_{\{t_3 \le t_4 \le t_1 \le t_2\}} \times (1-\lambda)^{t_2-t_3} N^{-2}. $$

We upper bound

$$ \sum_{t_1=1}^{\Delta} \sum_{t_2=1}^{\Delta} \sum_{t_3=1}^{\Delta} \sum_{t_4=1}^{\Delta} \mathbf{1}_{\{t_3 \le t_4 \le t_1 \le t_2\}} (1-\lambda)^{t_2-t_3} \le \Delta \left( \sum_{k=0}^{\Delta-1} (1-\lambda)^k \right)^3, $$

such that

$$ \sum_{\mathbf{z} \in C_{1,6}} \mathbf{1}_{\{t_3 \le t_4 \le t_1 \le t_2\}} |B(\mathbf{z})| \le K N^2 \Delta. $$

Proceeding similarly, we can also check that $\sum_{\mathbf{z} \in C_{1,6}} \mathbf{1}_{\{t_3 \le t_4 \le t_1 \le t_2\}} |B(\mathbf{z})| \le K N^2 \Delta$. This concludes the proof of Step 5.1

Step 5.2 Here we show that there exists a constant $K > 0$ such that for all $N \ge 1$ and $\Delta \ge 1$,

$$ \sum_{\mathbf{z} \in C_{1,6}} \left[ \mathbf{1}_{\{t_1 \le t_3 \le t_2 \le t_4\}} + \mathbf{1}_{\{t_3 \le t_1 \le t_4 \le t_2\}} + \mathbf{1}_{\{t_1 \le t_3 \le t_4 \le t_2\}} + \mathbf{1}_{\{t_3 \le t_1 \le t_2 \le t_4\}} \right] |B(\mathbf{z})| \le K N^2 \Delta^2. $$

From Item 7 of Lemma 3.8, we deduce that

$$ \sum_{\mathbf{z} \in C_{1,6}} \mathbf{1}_{\{t_1 \le t_3 \le t_2 \le t_4\}} |B(\mathbf{z})| \le K \sum_{\mathbf{z} \in C_{1,6}} \mathbf{1}_{\{t_1 \le t_3 \le t_2 \le t_4\}} (1-\lambda)^{(t_4-t_2)+(t_3-t_1)} N^{-2}. $$

Since

$$ \sum_{t_1=1}^{\Delta} \sum_{t_2=1}^{\Delta} \sum_{t_3=1}^{\Delta} \sum_{t_4=1}^{\Delta} \mathbf{1}_{\{t_1 \le t_3 \le t_2 \le t_4\}} (1-\lambda)^{(t_4-t_2)+(t_3-t_1)} \le \Delta^2 \left( \sum_{k=0}^{\Delta-1} (1-\lambda)^k \right)^2, $$

it then follows that

$$ \sum_{\mathbf{z} \in C_{1,6}} \mathbf{1}_{\{t_1 \le t_3 \le t_2 \le t_4\}} (1-\lambda)^{(t_1-t_3)+(t_2-t_4)} N^{-2} \le N^2 \Delta^2 \lambda^{-2} $$
such that

$$ \sum_{\mathbf{z} \in C_{1,6}} \mathbf{1}_{\{t_1 \le t_3 \le t_2 \le t_4\}} |B(\mathbf{z})| \le KN^2 \Delta^2. $$

Proceeding similarly, we can also check the three other cases, concluding the proof of Step 5.2.

Combining Steps 5.1, 5.2 and 5.3, we obtain

$$ \sum_{\mathbf{z} \in C_{1,6}} |B(\mathbf{z})| \le KN^2 \Delta^2, $$

concluding the proof of Step 5.

By summing all the upper bounds provided by the Steps 1 through 5, we conclude the proof of Lemma E.9. □

We are now in position to handle the terms $D_{\Delta,T}^{N,2}$ and $D_{2\Delta,T}^{N,2}$.

**PROPOSITION E.10.** There exists a constant $K$ depending only on $\lambda$ such that for all $N \ge 1$, $T \ge 2$ and $1 \le \Delta \le \lfloor T/2 \rfloor$, we have that

$$ \mathbb{E}_{\theta} \left[ \left( D_{\Delta,T}^{N,2} \right)^2 + \left( D_{2\Delta,T}^{N,2} \right)^2 \right] \le K \frac{\Delta}{T}. $$

**PROOF.** Clearly, the result will follow if we can show that $\mathbb{E}_{\theta} \left[ (D_{\Delta,T}^{N,2})^2 \right] \le K \Delta T^{-1}$. To establish this inequality, first observe that

$$ \begin{aligned} \mathbb{E}_{\theta} \left[ \left( D_{\Delta,T}^{N,2} \right)^2 \right] &= \frac{N^2}{T^2} \mathrm{Var}_{\theta} \left[ \sum_{t=1}^{\lfloor T/\Delta \rfloor} \left( \bar{U}_{t\Delta} - \bar{U}_{(t-1)\Delta} \right)^2 \right] \\ &= \frac{N^2}{T^2} \sum_{t=1}^{\lfloor T/\Delta \rfloor} \sum_{s=1}^{\lfloor T/\Delta \rfloor} \mathrm{Cov}_{\theta} \left[ \left( \bar{U}_{t\Delta} - \bar{U}_{(t-1)\Delta} \right)^2, \left( \bar{U}_{s\Delta} - \bar{U}_{(s-1)\Delta} \right)^2 \right]. \end{aligned} $$

Next, we use the stationarity of the system to conclude that

$$ \begin{aligned} \mathbb{E}_{\theta} \left[ \left( D_{\Delta,T}^{N,2} \right)^2 \right] &= \frac{N^2}{T^2} \lfloor T/\Delta \rfloor \mathrm{Var}_{\theta} \left[ (\bar{U}_{\Delta})^2 \right] \\ &\quad + 2 \frac{N^2}{T^2} \sum_{t=2}^{\lfloor T/\Delta \rfloor} \left( \lfloor T/\Delta \rfloor - t + 1 \right) \mathrm{Cov}_{\theta} \left[ (\bar{U}_{\Delta})^2, (\bar{U}_{t\Delta} - \bar{U}_{(t-1)\Delta})^2 \right]. \end{aligned} $$

As a consequence of Lemma E.8 and of the fact that $\sum_{k=0}^{\infty} (1-\lambda)^{k\Delta} < \infty$, it then follows that

$$ \mathbb{E}_{\theta} \left[ \left( D_{\Delta,T}^{N,2} \right)^2 \right] \le \frac{N^2}{T\Delta} \mathrm{Var}_{\theta} \left[ (\bar{U}_{\Delta})^2 \right] + K(\Delta T)^{-1}. $$

Hence, using Lemma E.9, we establish the desired inequality and conclude the proof. □
**E.3. Proof of Equation (33).** The proof relies on a decomposition of the error $|\hat{v} - v_{\infty}^N|$ (see Inequality (69)). Each term of the decomposition is then controlled by one of the four Propositions E.11-E.14. Their proofs are quite similar and less involved than the proofs given in the previous subsection.

Recall that

$$ \hat{v} = \frac{1}{T^2} \sum_{i=1}^{N} (Z_{i,T} - \bar{Z}_T)^2 - \frac{N}{T^2} \bar{Z}_T + \frac{1}{T^3} \sum_{i=1}^{N} (Z_{i,T})^2, $$

and that $\hat{m} = \bar{Z}_T/T$. Therefore,

$$ (69) \qquad |\hat{v} - v_{\infty}^N| \le \Upsilon_T^{N,1} + \Upsilon_T^{N,2} + \Upsilon_T^{N,3}, $$

where

$$ \Upsilon_T^{N,1} = \left| \sum_{i=1}^{N} \left[ \left( \frac{Z_{i,T}}{T} - \frac{\bar{Z}_T}{T} \right)^2 - (m_i^N - \overline{m^N})^2 \right] - \frac{N}{T} \overline{m^N} + \frac{1}{T} \sum_{i=1}^{N} (m_i^N)^2 \right|, $$

$$ \Upsilon_T^{N,2} = \frac{N}{T} |\hat{m} - \overline{m^N}|, \quad \text{and} \quad \Upsilon_T^{N,3} = \frac{1}{T} \left| \sum_{i=1}^{N} \left( \frac{Z_{i,T}}{T} \right)^2 - (m_i^N)^2 \right|. $$

**PROPOSITION E.11.** There exists a constant K depending only on $\lambda$ such that for all $N \ge 1$ and $T \ge 1$,

$$ \mathbb{E}_{\theta} \left[ \Upsilon_{T}^{N,2} + \Upsilon_{T}^{N,3} \right] \le K \left( \frac{N}{T^2} + \frac{N^{1/2}}{T^{3/2}} \left( 1 + \sqrt{v_{\infty}^N} \right) \right). $$

**PROOF.** First, remind that $\mathbb{E}_{\theta}[Z_{i,T}/T] = m_i^N$ and so $\mathbb{E}_{\theta}[\hat{m}] = \overline{m^N}$. On the one hand, by Cauchy-Schwarz, $\mathbb{E}_{\theta}[\Upsilon_T^{N,2}] \le NT^{-1}(\mathrm{Var}_{\theta}(\hat{m}))^{1/2} \le KN^{1/2}T^{-3/2}$ by Equation (32). On the other hand, by first rewriting $\Upsilon_T^{N,3}$ as (recall that $U_{i,T} = Z_{i,T} - Tm_i^N$),

$$ \Upsilon_T^{N,3} = \frac{1}{T} \left| \sum_{i=1}^{N} \left( \frac{U_{i,T}}{T} \right)^2 + 2 \sum_{i=1}^{N} (m_i^N - \overline{m^N}) \frac{U_{i,T}}{T} + 2 \frac{N}{T} \overline{m^N} \bar{U}_T \right| $$

and then using the triangle inequality, we obtain that

$$ \mathbb{E}_{\theta} \left[ \Upsilon_{T}^{N,3} \right] \le \frac{1}{T^3} \sum_{i=1}^{N} \mathbb{E}_{\theta} \left[ (U_{i,T})^2 \right] + \frac{2}{T^2} \sum_{i=1}^{N} |m_i^N - \overline{m^N}| \mathbb{E}_{\theta} \left[ |U_{i,T}| \right] + 2 \frac{N}{T^2} \overline{m^N} \mathbb{E}_{\theta} \left[ |\bar{U}_T| \right]. $$

Now, by stationarity and Lemma 3.7,

$$ (70) \qquad \begin{aligned} \mathrm{Var}_{\theta}(Z_{i,T}) &= T \mathrm{Var}_{\theta}(X_{i,0}) + 2 \sum_{t=1}^{T} (T-t) \mathrm{Cov}_{\theta}(X_{i,0}, X_{i,t}) \\ &\le T \mathrm{Var}_{\theta}(X_{i,0}) + KT/N \le KT, \end{aligned} $$

so that

$$ \frac{1}{T^3} \sum_{i=1}^{N} \mathbb{E}_{\theta} \left[ (U_{i,T})^2 \right] = \frac{1}{T^3} \sum_{i=1}^{N} \mathrm{Var}_{\theta}(Z_{i,T}) \le KNT^{-2}. $$
Moreover, by using Cauchy-Schwarz inequality and then Jensen inequality, we deduce that

$$ \sum_{i=1}^{N} |m_{i}^{N} - \overline{m}^{N}| \mathcal{E}_{\theta} [|U_{i,T}|] \leq \sqrt{\sum_{i=1}^{N} |m_{i}^{N} - \overline{m}^{N}|^{2}} \sqrt{\sum_{i=1}^{N} \mathcal{E}_{\theta} [(U_{i,T})^{2}]} $$

Hence, using once more Inequality (70) and reminding that $v_{\infty}^{N} = \sum_{i=1}^{N} |m_{i}^{N} - \overline{m}^{N}|^{2}$, we obtain that

$$ \frac{2}{T^2} \sum_{i=1}^{N} |m_i^N - \overline{m}^N| \mathcal{E}_{\theta} [|U_{i,T}|] \leq K \frac{N^{1/2}}{T^{3/2}} \sqrt{v_{\infty}^N} $$

Finally, by combining Jensen inequality and Lemma E.6, it follows that $\mathcal{E}_{\theta} [|\overline{U}_{T}|] \leq \sqrt{\mathcal{E}_{\theta} [|\overline{U}_{T}|^2]} \leq KT^{1/2}N^{-1/2}$ so that

$$ 2 \frac{N}{T^2} \overline{m}^N \mathcal{E}_{\theta} [|\overline{U}_T|] \leq K N^{1/2} T^{-3/2}, $$

where we have also used that $|\overline{m}^N| \leq 1$. Putting together the above estimates, it then follows that

$$ \mathcal{E}_{\theta}[\Upsilon_{T}^{N,3}] \leq K \left( \frac{N}{T^2} + \frac{N^{1/2}}{T^{3/2}} \sqrt{v_{\infty}^N} \right), $$

and the result follows since we have already proved that $\mathcal{E}_{\theta}[\Upsilon_{T}^{N,2}] \leq KN^{1/2}T^{-3/2}$. $\square$

<figure><img src="image_15.png" /></figure>

We now turn to the study of $\Upsilon_{T}^{N,1}$. Using that

$$ (a-b)^2 - (\bar{a} - \bar{b})^2 = (a-b)^2 - (a-\bar{b})^2 + (a-\bar{a})^2 + 2(a-\bar{a})(\bar{a}-\bar{b}), $$

we obtain that $\Upsilon_{T}^{N,1} \leq \Upsilon_{T}^{N,1,1} + \Upsilon_{T}^{N,1,2} + \Upsilon_{T}^{N,1,3}$, where

$$ \Upsilon_{T}^{N,1,1} = \left| \sum_{i=1}^{N} \left( \frac{Z_{i,T}}{T} - \frac{\bar{Z}_{T}}{T} \right)^{2} - \left( \frac{Z_{i,T}}{T} - \overline{m}^{N} \right)^{2} \right|, $$

$$ \Upsilon_{T}^{N,1,2} = \left| \sum_{i=1}^{N} \left( \frac{Z_{i,T}}{T} - m_{i}^{N} \right)^{2} - \frac{N}{T} \overline{m}^{N} + \frac{1}{T} \sum_{i=1}^{N} (m_{i}^{N})^{2} \right| $$

and

$$ \Upsilon_{T}^{N,1,3} = 2 \left| \sum_{i=1}^{N} \left( \frac{Z_{i,T}}{T} - m_{i}^{N} \right) (m_{i}^{N} - \overline{m}^{N}) \right|. $$

Furthermore, remind that $U_{i,t} = Z_{i,t} - tm_{i}^{N}$ and so we write $\Upsilon_{T}^{N,1,2} = \Upsilon_{T}^{N,1,2,1} + \Upsilon_{T}^{N,1,2,2}$, with

$$ \Upsilon_{T}^{N,1,2,1} = \left| \sum_{i=1}^{N} \left( \frac{U_{i,T}}{T} \right)^{2} - \mathcal{E}_{\theta} \left( \left( \frac{U_{i,T}}{T} \right)^{2} \right) \right| $$

and

$$ \Upsilon_{T}^{N,1,2,2} = \left| \sum_{i=1}^{N} \mathcal{E}_{\theta} \left( \left( \frac{U_{i,T}}{T} \right)^{2} \right) - \frac{N}{T} \overline{m}^{N} + \frac{1}{T} \sum_{i=1}^{N} (m_{i}^{N})^{2} \right|. $$
PROPOSITION E.12. For all $N \ge 1$ and $T \ge 1$, there exists a constant $K$ depending only on $\lambda$ such that

$$ \mathbb{E}_{\theta} \left[ \Upsilon_{T}^{N,1,1} + \Upsilon_{T}^{N,1,2,2} \right] \le K \frac{1}{T}. $$

PROOF. Adapting the argument used in the beginning of Proposition E.1, we have

$$ \Upsilon_{T}^{N,1,1} = N \left( \frac{\bar{Z}_{T}}{T} - \overline{m^{N}} \right)^{2}. $$

Hence, Equation (32) implies that $\mathbb{E}_{\theta}\Upsilon_{T}^{N,1,1} \le K/T$.

Second, remark that $\Upsilon_{T}^{N,1,2,2}$ is not random and rewrites as

$$ \Upsilon_{T}^{N,1,2,2} = \left| \sum_{i=1}^{N} \mathrm{Var}_{\theta} \left( \frac{Z_{i,T}}{T} \right) - \frac{N}{T} \overline{m^{N}} + \frac{1}{T} \sum_{i=1}^{N} (m_{i}^{N})^{2} \right|. $$

Using Equation (70) and the fact that $\mathrm{Var}_{\theta}(X_{i,0}) = m_{i}^{N} - (m_{i}^{N})^{2}$ we conclude that $\mathbb{E}_{\theta}\Upsilon_{T}^{N,1,2,2} \le K/T$. $\square$

PROPOSITION E.13. For all $N \ge 1$ and $T \ge 1$, there exists a constant $K$ depending only on $\lambda$ such that

$$ \mathbb{E}_{\theta}|\Upsilon_{T}^{N,1,3}|^2 \le \frac{K}{T} \sum_{i=1}^{N} |m_i^N - \overline{m^N}|^2. $$

PROOF. Clearly,

$$ \mathbb{E}_{\theta}|\Upsilon_{T}^{N,1,3}|^2 \le \frac{4}{T^2} \sum_{i,j=1}^{N} |\mathrm{Cov}_{\theta}(U_{i,T}, U_{j,T})| |m_i^N - \overline{m^N}| |m_j^N - \overline{m^N}|. $$

But, by stationarity and using Lemma 3.7,

$$ |\mathrm{Cov}_{\theta}(U_{i,T}, U_{j,T})| \le \sum_{s,t=1}^{T} |\mathrm{Cov}_{\theta}(X_{i,s}, X_{j,t})| \le T |\mathrm{Cov}_{\theta}(X_{i,0}, X_{j,0})| + KT/N. $$

Moreover,

$$ |\mathrm{Cov}_{\theta}(X_{i,0}, X_{j,0})| \le (m_i^N - (m_i^N)^2) \mathbf{1}_{\{i=j\}} + K/N \mathbf{1}_{\{i \neq j\}}, $$

such that all in all

$$ \mathbb{E}_{\theta}|\Upsilon_{T}^{N,1,3}|^2 \le \frac{K}{T} \left( \sum_{i=1}^{N} \left[ \frac{1}{N} + (m_i^N - (m_i^N)^2) \right] (m_i^N - \overline{m^N})^2 + \frac{1}{N} \sum_{i \neq j} |m_i^N - \overline{m^N}| |m_j^N - \overline{m^N}| \right). $$

Since $m_i^N - (m_i^N)^2 \le 1$ and $|m_i^N - \overline{m^N}| |m_j^N - \overline{m^N}| \le (m_i^N - \overline{m^N})^2 + (m_j^N - \overline{m^N})^2$, the conclusion follows. $\square$

PROPOSITION E.14. For all $N \ge 1$ and $T \ge 1$, there exists a constant $K$ depending only on $\lambda$ such that

$$ \mathbb{E}_{\theta}|\Upsilon_{T}^{N,1,2,1}|^2 \le KN/T^2. $$
PROOF. We have that

$$
\begin{align*}
(71) \quad \mathbb{E}_{\theta} |\Upsilon_T^{N,1,2,1}|^2 &= \frac{1}{T^4} \sum_{i,j=1}^{N} \mathrm{Cov}_{\theta}((U_{i,T})^2, (U_{j,T})^2) \\
&= \frac{1}{T^4} \sum_{i,j=1}^{N} \sum_{s=1}^{T} \sum_{t=1}^{T} \mathrm{Cov}_{\theta}(Y_{i,s} Y_{i,s}, Y_{j,t} Y_{j,t}) + \frac{1}{T^4} \sum_{i,j=1}^{N} \sum_{s \neq s'} \sum_{t \neq t'} \mathrm{Cov}_{\theta}(Y_{i,s} Y_{i,s'}, Y_{j,t} Y_{j,t'}) \\
&=: S_1 + S_2.
\end{align*}
$$

By Lemma 3.8, Item 1 and Item 3, $S_1 \le KN/T^3$. To deal with $S_2$, we adapt the arguments of the proof of Lemma E.9. As before, write $\mathbf{z} = (z_1, z_2, z_3, z_4)$, and $B(\mathbf{z}) = \mathrm{Cov}_{\theta}[Y_{z_1} Y_{z_2}, Y_{z_3} Y_{z_4}]$. This expression is symmetric with respect to $t, t'$ and also with respect to $s, s'$. Now, introduce the set $F^{(2)} = \{(z_1, z_2) : z_1 = (i, s), z_2 = (i, s'), i \in \{1, \dots, N\}, s, s' \in \{1, \dots, T\}, s \neq s'\}$ and define, for $p \in \{2, 5, 7\}$, $P_{p,T}$ as the set of all vectors $\mathbf{z} \in F^{(2)} \times F^{(2)}$ satisfying the conditions of Item $p$ of Lemma 3.8. Moreover, notice that $\mathbf{z} \in P_{5,T} \cup P_{2,T}$ implies that the first coordinates (denoted $i$ and $j$ above) of the four couples in $\mathbf{z}$ are equal.

Then

$$S_2 = \frac{1}{T^4} \left[ \sum_{\mathbf{z} \in \mathcal{P}_{7,T}} B(\mathbf{z}) + \sum_{\mathbf{z} \in \mathcal{P}_{5,T}} B(\mathbf{z}) + \sum_{\mathbf{z} \in \mathcal{P}_{2,T}} B(\mathbf{z}) \right].$$

Following Step 5 of the proof of Lemma E.9, it is easy so see that

$$\sum_{\mathbf{z} \in \mathcal{P}_{7,T}} B(\mathbf{z}) \le KT^2.$$

Moreover, Lemma 3.8, Item 5, implies that

$$\sum_{\mathbf{z} \in \mathcal{P}_{5,T}} B(\mathbf{z}) \le KT^2.$$

Finally,

$$\sum_{\mathbf{z} \in \mathcal{P}_{2,T}} B(\mathbf{z}) \le K |\mathcal{P}_{2,T}| = KNT^2.$$

All in all we therefore obtain that $S_2 \le KN/T^2$, implying the assertion. $\square$

## APPENDIX F: PROOF OF PROPOSITION 4.2

This section contains a fine study of the random environment $\theta$. First, some matrix notation is introduced. Then, the asymptotics of the rows and columns of the rescaled random environment $A^N$ are stated in Lemma F.1. In turn, the asymptotics of the rows and columns of the inverse matrix $Q^N$ are stated in Lemmas F.4 and F.6. Then, the proofs of these three lemmas are given in three disjoint subsections. Finally, the proof of Proposition 4.2 is given in the last subsection.

To readers interested in random environments that differ from the i.i.d. case, let us mention that the proof of Proposition 4.2 relies solely on the results of Lemmas F.1, F.4, and F.6 in the sense that the same proof applies to any environment $\theta$ for which these lemmas hold. Similarly, the proofs of Lemmas F.4 and F.6 depend only on Lemma F.1 Therefore, only the proof Lemma F.1 must be adapted. In what follows, we specify the adaptation required for the symmetric interaction case.
The main modification concerns the bound of the variance terms for which we use Bien-aymé's identity together with a control of the few non null covariance terms. The exhaustive list of modifications is: 1) Equation (96), 2) Equation (97) by proving that $\mathrm{Var}(\sum_{k \in \mathcal{P}_a} Z_k)$ and $\mathrm{IE}[(\sum_{k \in \mathcal{P}_a} \theta_{1j} Z_j)^2]$ are smaller than some universal constant, 3) check that the "independence argument" below Equation (98) still holds, 4) replace the "independence argument" by Cauchy-Schwarz inequality to bound $\mu_{1,1}$, 5) Equation (99) is replaced by

$$ \mathrm{IE} \left[ \left( Z_i(1 - |\mathcal{P}_a|^{-1}) - \frac{1}{|\mathcal{P}_a|} \sum_{j \in \mathcal{P}_a: j \neq i} Z_j \right)^2 \right] \leq (1 - |\mathcal{P}_a|^{-1}) \mathrm{Var}(Z_i) + KN^{-2} $$

for some universal constant $K$.

Two technical results used throughout the proofs are stated as technical lemmas in Section H.

**F.1. General notation.** Hereafter, for any subset $S$ of $[N]$, we write $1_S$ to indicate the $N$-dimensional vector having value 1 in each coordinate belonging to $S$ and value 0 in the remaining coordinates. To alleviate the notation, we will simply write $1_N$ and $0_N$ instead of $1_{[N]}$ and $1_\emptyset$, respectively. For any vector $v \in \mathbb{R}^N$, $\bar{v} = N^{-1} \sum_{i=1}^N v_i$ denotes the arithmetic mean of the coordinates of $v$, $\|v\|_r = (\sum_{i=1}^N |v_i|^r)^{1/r}$ where $r \in [1, \infty)$, denotes the $r$-norm of $v$ and $\|v\|_\infty = \max_{1 \leq i \leq N} |v_i|$ its $\infty$-norm. For vectors $v, u \in \mathbb{R}^N$, we write $v \odot u$ to denote the vector whose $i$-th coordinate is $v_i u_i$, for $1 \leq i \leq N$. In other words, $v \odot u$ is the Hadamard product between the vectors $v$ and $u$. To shorten the notation, we will simply note $v^2$ instead of $v \odot v$. Observe that, with this notation, the orthogonal projection of a vector $v \in \mathbb{R}^N$ on the coordinates in $S \subseteq [N]$ can be written as $1_S \odot v$. In particular, all coordinates in $S^c$ of the vector $1_S \odot v$ are null. One can always write $v = 1_S \odot v + 1_{S^c} \odot v$.

For any $N$-by-$N$ matrix $B$ with real entries, we denote $B^\top$ its transpose. For all $r \in [1, \infty]$, we denote $\|B\|_r$ the operator norm of $B$ associated to the $r$-norm $\| \cdot \|_r$:

$$ \|B\|_r = \sup_{v \in \mathbb{R}^N: v \neq 0_N} \frac{\|Bv\|_r}{\|v\|_r}. $$

It is well-known that $\|B\|_1$ and $\|B\|_\infty$ may be defined alternatively as

$$ \|B\|_1 = \max_{1 \leq j \leq N} \sum_{i=1}^N |B(i, j)| \quad \text{and} \quad \|B\|_\infty = \max_{1 \leq i \leq N} \sum_{j=1}^N |B(i, j)|. $$

The following fact will also be used in the sequel: for each $r \in (1, \infty)$, it holds that

$$ (72) \qquad \|B\|_r \leq \|B\|_1^{1/r} \|B\|_\infty^{1-1/r}. $$

**F.2. Study of the rescaled random environment.** Recall that $A^N = (A^N(i, j))_{1 \leq i, j \leq N}$ is a rescaled version of the random environment $\theta$ defined as, for each $1 \leq i \leq N$,

$$ A^N(i, j) = \begin{cases} N^{-1}\theta_{ij}, & \text{if } j \in \mathcal{P}_+, \\ -N^{-1}\theta_{ij}, & \text{if } j \in \mathcal{P}_-. \end{cases} $$

One can check that $\max\{\|A^N\|_1, \|A^N\|_\infty\} \leq 1$, so that (72) implies that $\|A^N\|_r \leq 1$ for all $r \in (1, \infty)$ as well.
F.2.1. Notation. In what follows, for $a, b \in \{-, +\}$, we denote $L^{N,ab} = 1_{P_a} \odot A^N 1_{P_b}$ and $C^{N,ab} = 1_{P_a} \odot (A^N)^\top 1_{P_b}$, where $(A^N)^\top$ denotes the transpose of the matrix $A^N$. On the one hand, both vectors have null coordinates outside of the set $P_a$. On the other hand, each coordinate $i \in P_a$ of the random vector $L^{N,ab}$ (resp. $C^{N,ab}$) is obtained by summing the entries in $P_b$ of the $i$-th row (resp. column) of the matrix $A^N$. Alternatively, the random vectors $L^{N,ab}$ and $C^{N,ab}$ can be defined as follows:

$$ (73) \qquad L_i^{N,ab} = \begin{cases} \sum_{j \in P_b} A^N(i, j), & \text{if } i \in P_a, \\ 0, & \text{otherwise,} \end{cases} $$

and

$$ (74) \qquad C_i^{N,ab} = \begin{cases} \sum_{j \in P_b} A^N(j, i), & \text{if } i \in P_a, \\ 0, & \text{otherwise.} \end{cases} $$

For $a \in \{-, +\}$, we denote $L^{N,a\bullet} = 1_{P_a} \odot A^N 1_N$ and $C^{N,a\bullet} = 1_{P_a} \odot (A^N)^\top 1_N$. Note that the coordinates not belonging to $P_a$ of these two vectors are also 0. Besides, each coordinate $i \in P_a$ of the random vector $L^{N,a\bullet}$ (resp. $C^{N,a\bullet}$) is given by the sum over all entries of the $i$-th row (resp. column) of the matrix $A^N$. One can easily check that

$$ (75) \qquad L^{N,a\bullet} = L^{N,a+} + L^{N,a-} \quad \text{and} \quad C^{N,a\bullet} = C^{N,a+} + C^{N,a-}, $$

so that the vectors $L^{N,a\bullet}$ and $C^{N,a\bullet}$ could be defined alternatively through these identities.

For $b \in \{-, +\}$, we denote $L^{N,b\bullet} = A^N 1_{P_b}$ and $C^{N,b\bullet} = (A^N)^\top 1_{P_b}$. Observe that each coordinate $1 \le i \le N$ of the random vector $L^{N,b\bullet}$ (resp. $C^{N,b\bullet}$) is given by the sum over the entries in $P_b$ of the $i$-th row (resp. column) of the matrix $A^N$. One can also verify that

$$ (76) \qquad L^{N,b\bullet} = L^{N,+b} + L^{N,-b} \quad \text{and} \quad C^{N,b\bullet} = C^{N,+b} + C^{N,-b}. $$

In the sequel, for each $a, b \in \{-, +\}$, let

$$ (77) \qquad \tilde{L}^{N,ab} = L^{N,ab} - (r_a^N)^{-1} \overline{L^{N,ab}} 1_{P_a} \quad \text{and} \quad \tilde{C}^{N,ab} = C^{N,ab} - (r_a^N)^{-1} \overline{C^{N,ab}} 1_{P_a}. $$

The vectors $\tilde{L}^{N,ab}$ and $\tilde{C}^{N,ab}$ can be thought as the population-wise centered versions of the vectors $L^{N,ab}$ and $C^{N,ab}$ respectively, in the sense that $\tilde{L}^{N,ab} = \tilde{C}^{N,ab} = 0$. Observe that these vectors are well-defined for all $N \ge 1$ such that $r_a^N > 0$. In the next result, we collect some bounds which will be used throughout the section.

F.2.2. Convergence rates. In what follows, for $a, b \in \{-, +\}$, we write $\delta_{ab}$ to denote the Kronecker delta between $a$ and $b$.

LEMMA F.1. Let $a, b, a_1, b_1 \in \{-, +\}$. The following inequalities hold for all $N \ge 1$ such that $r_{-}^N \wedge r_{+}^N > 0$,

$$ (78) \qquad \mathrm{IE} \left[ \left| \langle \tilde{L}^{N,ab}, \tilde{L}^{N,a_1 b_1} \rangle - \delta_{a_1 a} \delta_{b_1 b} r_a r_b p(1-p) \right|^2 \right] \le KN^{-1}, $$

$$ (79) \qquad \mathrm{IE} \left[ \| A \tilde{L}^{ab} \|_2^2 \right] \le KN^{-1}, $$

$$ (80) \qquad \mathrm{IE} \left[ \max_{a,b \in \{-, +\}} \left\{ \left| (r_a^N r_b^N)^{-1} \overline{L^{N,ab}} - (bp)^2 \right| \vee \left| (r_a^N r_b^N)^{-1} \overline{C^{N,ab}} - (ap)^2 \right| \right\} \right] \le K \max_{a,b} \{(r_a^N r_b^N)^{-1}\} N^{-2}, \quad \text{and} $$

$$ (81) \qquad \| \tilde{L}^{N,ab} \|_2 \le \| L^{N,ab} - (bpr_b^N) 1_{P_a} \|_2 \quad \text{and} \quad \| \tilde{C}^{N,ab} \|_2 \le \| C^{N,ab} - (apr_b^N) 1_{P_a} \|_2, $$
where $K$ is some universal constant. Moreover, for each $\alpha \ge 1$, there exists a constant $K_\alpha$ depending only on $\alpha$ such that, for all $N \ge 1$ satisfying $r_-^N \wedge r_+^N > 0$,

$$ (82) \quad \mathbb{E} \left[ \max_{a,b \in \{-, +\}} \left\{ \|L^{N,ab} - (bpr_b^N)1_{\mathcal{P}_a}\|_2^\alpha \vee \|C^{N,ab} - (apr_b^N)1_{\mathcal{P}_a}\|_2^\alpha \right\} \right] \\ \le K_\alpha \max_{a \in \{-, +\}} \left\{ (r_a^N/r_b^N)^{\alpha/2} \right\}. $$

The proof of Lemma F.1 is postponed to Section F.4.

REMARK F.2. Let $a, b \in \{-, +\}$. Because of the relations (75) and (76), some statements of Lemma F.1 admit immediate corollaries. For instance,

* $C^{N, \bullet b} = C^{N,+b} + C^{N,-b}$ and $\tilde{C}^{N, \bullet b} = \tilde{C}^{N,+b} + \tilde{C}^{N,-b}$ are orthogonal sums, and in particular one can deduce from Equations (81) and (82) that

$$ \mathbb{E} \left[ \| \tilde{C}^{N, \bullet b} \|_2^\alpha \right] \le \mathbb{E} \left[ \| C^{N, \bullet b} - pr_b^N (1_{\mathcal{P}_+} - 1_{\mathcal{P}_-}) \|_2^\alpha \right] \le 2K_\alpha \max_{a \in \{-, +\}} \left\{ (r_a^N/r_b^N)^{\alpha/2} \right\}, $$

* or, $L^{N,a\bullet} = L^{N,a+} + L^{N,a-}$ and $\tilde{L}^{N,a\bullet} = \tilde{L}^{N,a+} + \tilde{L}^{N,a-}$, and in particular one can deduce from Equations (81) and (82) that

$$ \|\tilde{L}^{N,a\bullet}\|_2 \le \|L^{N,a+} - pr_+^N 1_{\mathcal{P}_a}\|_2 + \|L^{N,a-} + pr_-^N 1_{\mathcal{P}_a}\|_2, $$

and

$$ \mathbb{E} \left[ \| L^{N,a\bullet} - p(r_+^N - r_-^N) 1_{\mathcal{P}_a} \|_2^\alpha \right] \le K'_\alpha \max_{a \in \{-, +\}} \left\{ (r_a^N/r_b^N)^{\alpha/2} \right\}, $$

for some $K'_\alpha$ which may be different from $K_\alpha$.

F.3. **Results regarding the inverse matrix.** Recall that $\lambda > 0$. Under this condition, the random matrix $Q^N := (I - (1 - \lambda)A^N)^{-1} = \sum_{n=0}^\infty (1 - \lambda)^n (A^N)^n$ is well-defined and satisfies $\|Q^N\|_r \le \lambda^{-1}$ for all $r \in [1, \infty]$. For later use, let us observe that the easy-to-check properties $\|(Q^N)^\top\|_1 = \|Q^N\|_\infty$ and $\|(Q^N)^\top\|_\infty = \|Q^N\|_1$ combined with inequality (72) imply that $\|(Q^N)^\top\|_r \le \lambda^{-1}$ for any $r \in [1, \infty]$.

As suggested by the heuristics presented in Section 2, a crucial ingredient in our analysis is to study both the sum of rows and columns of the matrix $Q^N$. By definition of the matrix $Q^N$ itself, these quantities are related to the corresponding counterparts computed from the matrix $A^N$.

F.3.1. **Notation.** Similarly as above, for $a, b \in \{-, +\}$, let us denote $\ell^{N,ab}$ (resp. $c^{N,ab}$) the random vector obtained by summing, for each row (resp. column) in $\mathcal{P}_a$ of the random matrix $Q^N$, the entries in $\mathcal{P}_b$. Also, we denote $\ell^{N,a\bullet} = \ell^{N,a+} + \ell^{N,a-}$ and $\tilde{\ell}^{N,\bullet b} = \ell^{N,+b} + \ell^{N,-b}$ for $a, b \in \{-, +\}$, and define $\ell^N = \ell^{N,+\bullet} + \ell^{N,-\bullet}$. Note that $\ell^N = Q^N 1_N$, i.e., $\ell^N$ corresponds to the random vector obtained by summing the rows of $Q^N$. The vectors $c^{N,\bullet a}$, $c^{N,\bullet a}$ and $c^N$ are defined in a similar way. In particular, note that $c^N = (Q^N)^\top 1_N$ is given as the sum of the columns of $Q^N$. Here $(Q^N)^\top$ denotes the transpose of the matrix $Q^N$. For later use, let us also observe that $\overline{\ell^N} = \overline{c^N}$.

Like we did for $L$ and $C$, we define, for each $a, b \in \{-, +\}$

$$ (83) \quad \tilde{\ell}^{N,ab} = \ell^{N,ab} - (r_a^N)^{-1} \overline{\ell^{N,ab}} 1_{\mathcal{P}_a} \quad \text{and} \quad \tilde{c}^{N,ab} = c^{N,ab} - (r_a^N)^{-1} \overline{c^{N,ab}} 1_{\mathcal{P}_a}. $$

Moreover, we denote $\tilde{\ell}^{N,a\bullet} = \tilde{\ell}^{N,a+} + \tilde{\ell}^{N,a-}$, $\tilde{\ell}^{N,\bullet b} = \tilde{\ell}^{N,+b} + \tilde{\ell}^{N,-b}$ for $a, b \in \{-, +\}$, and $\overline{\ell^N} = \overline{\ell^{N,+\bullet} + \ell^{N,-\bullet}}$.

Here is an immediate result stating that the $\ell$ and $c$ vectors are uniformly bounded.
LEMMA F.3. Assume that $1 \ge \lambda > 0$. Then,

$$ (84) \quad \max \left\{ \max_{a,b \in \{-, +\}} \left\{ \|\ell^{N,ab}\|_{\infty}, \|c^{N,ab}\|_{\infty} \right\}, \|\ell^N\|_{\infty}, \|c^N\|_{\infty} \right\} \le \lambda^{-1}. $$

PROOF. To show (84), first observe that

$$ \begin{gather*} \max_{1 \le i \le N} \max_{a,b \in \{-, +\}} \left\{ |\ell_i^{N,ab}|, |c_i^{N,ab}| \right\}, |\ell_i^N|, |c_i^N| \right\} \le \max\{\|Q^N\|_1, \|Q^N\|_{\infty}\}, \\ \text{and then use the fact that } \max\{\|Q^N\|_1, \|Q^N\|_{\infty}\} \le \lambda^{-1}. \quad \square \end{gather*} $$

F.3.2. *Convergence rates.* Recall that $\delta_{ab}$ denotes the Kronecker delta between $a$ and $b$.

LEMMA F.4. Assume that $1 \ge \lambda > 0$. There exists a constant $K > 0$ which depends on $\lambda$ such that for all $N \ge 1$ and $a, b \in \{-, +\}$, it holds that

$$ (85) \quad \mathbb{E} \left[ \left| \overline{\ell^{N,ab}} - \frac{r_a \delta_{ab} + b(1-\lambda)pr_+r_-}{1-(1-\lambda)p(r_+-r_-)} \right|^2 \right] \le \frac{K}{N^2}, $$

$$ (86) \quad \mathbb{E} \left[ \left| c^{N,a\bullet} - r_a \frac{[1+(2a)(1-\lambda)p(1-r_a)]}{1-(1-\lambda)p(r_+-r_-)} \right|^2 \right] \le \frac{K}{N^2}, $$

$$ (87) \quad \mathbb{E} \left[ \| \tilde{\ell}^{N, \bullet b} \|_2^4 \right] \le K \quad \text{and} \quad \mathbb{E} \left[ \| \tilde{c}^{N, \bullet b} \|_2^4 \right] \le K, $$

$$ (88) \quad \mathbb{E} \left[ \left\| \tilde{\ell}^{N, \bullet b} \right\|_2^2 - \frac{(1-\lambda)^2 p(1-p) [r_b + (1-\lambda)pr_+r_-(2b+(1-\lambda)p)]}{(1-(1-\lambda)p(r_+-r_-))^2} \right] \le \frac{K}{\sqrt{N}}, $$

$$ (89) \quad \mathbb{E} \left[ \left\| \tilde{\ell}^N \right\|_2^2 - \frac{(1-\lambda)^2 p(1-p)}{(1-(1-\lambda)p(r_+-r_-))^2} \right] \le \frac{K}{\sqrt{N}}. $$

REMARK F.5. First, notice that Remark F.2 also applies to Lemma F.4 (with $L$ and $C$ replaced by $\ell$ and $c$).

For later use, let us observe that an immediate consequence of Equation (85) is that there exists a constant $K > 0$ such that $N \ge 1$,

$$ (90) \quad \mathbb{E} \left[ \left| \overline{\ell^N} - \frac{1}{1-(1-\lambda)p(r_+-r_-)} \right|^2 \right] \le \frac{K}{N^2}. $$

The proof of Lemma F.4 is postponed to Section F.5.

The following lemma gives a way to control the $l^2$ norms of the fully centered vectors by the population wise centered vectors.

LEMMA F.6. Assume that $1 \ge \lambda > 0$. There exists a constant $K > 0$ which depends on $\lambda$ such that for all $N \ge 1$ such that $r_+^N \wedge r_-^N > 0$, it holds that

$$ (91) \quad \mathbb{E} \left[ \left\| \ell^N - \overline{\ell^N} 1_N \right\|_2^2 - \left\| \tilde{\ell}^N \right\|_2^2 \right] \le \frac{K}{N}, $$

$$ (92) \quad \mathbb{E} \left[ \left| \langle \ell^N - \overline{\ell^N} 1_N, (1-\lambda)^2 (Q^N L^{N,\bullet} - \overline{Q^N L^{N,\bullet}}) \rangle - \langle \tilde{\ell}^N, \tilde{\overline{\ell^N}} \rangle \right| \right] \le KN^{-1}, $$

$$ (93) \quad \mathbb{E} \left[ \left| (1-\lambda) \|QL^{N,\bullet} - \overline{Q^N L^{N,\bullet}}\|_2^2 - \|\tilde{\ell}^{N,\bullet} - \overline{\tilde{\ell^N}}\|_2^2 \right| \right] \le KN^{-1}. $$

The proof of Lemma F.6 is postponed to Section F.6.
**F.4.** *Proof of Lemma F.1.* Before the proof, we recall some classical results in the three lemmas below.

**LEMMA F.7.** Consider i.i.d. random variables $(B_i)_{1\le i\le N}$ distributed as Ber($p$) with $p \in [0, 1]$ and let $\bar{B} = N^{-1}\sum_{i=1}^N B_i$. Then, for each $\alpha \ge 1$ there exists a constant $C_\alpha > 0$ depending only on $\alpha$ such that

$$ \mathbb{E}[|\bar{B} - p|^\alpha] \le C_\alpha N^{-\alpha/2}. $$

**PROOF.** Hoeffding's inequality implies that for any $x > 0$,

$$ \mathbb{P}(|\bar{B} - p| \ge x) \le 2e^{-2x^2N}. $$

Now, using the tail sum formula for the expectation one can show that

$$ \mathbb{E}[|\bar{B} - p|^\alpha] \le 2\alpha \int_0^\infty x^{\alpha-1}e^{-2x^2N}dx = \alpha\sqrt{\frac{2\pi}{N}}\mathbb{E}[Z^{\alpha-1}1_{Z>0}], $$

where $Z \sim \mathcal{N}(0, (4N)^{-1})$. Finally, by observing that $Y = (4N)^{1/2}Z \sim \mathcal{N}(0, 1)$, it follows that

$$ \alpha\sqrt{\frac{2\pi}{N}}\mathbb{E}[Z^{\alpha-1}1_{Z>0}] = \alpha\frac{\sqrt{2\pi}}{2^{\alpha-1}}N^{-\alpha/2}\mathbb{E}[Y^{\alpha-1}1_{Y>0}], $$

and the result follows. $\square$

**LEMMA F.8.** Let $(B_{ij})_{1\le i,j\le N}$ be i.i.d. random variables distributed as Ber($p$) with $p \in [0, 1]$. For $U, V \subset [N]$ and $\epsilon \in \{-1, 1\}$, define $W_i = \epsilon N^{-1} \sum_{j \in V} B_{ij}$ for $i \in U$, and $W_i = 0$, for $i \in U^c$. Denote $W = (W_1, \dots, W_N)$. Then, for any $\alpha \ge 1$, there exists a constant $K_\alpha > 0$ depending only on $\alpha$ such that

$$ (94) \qquad \mathbb{E} \left( \|W - p\epsilon |V| N^{-1} 1_U\|_2^\alpha \right) \le K_\alpha \left( \frac{|U|}{|V|} \right)^{\alpha/2}. $$

**PROOF.** Jensen inequality combined with Lemma F.7 and the fact that $|V|/N \le 1$ implies that for any $\alpha \ge 1$,

$$ \mathbb{E} \left[ \|W - p\epsilon |V| N^{-1} 1_U \|_2^{2\alpha} \right] \le |U|^{\alpha} (|V| N^{-1})^{2\alpha} (C_{2\alpha} |V|^{-(2\alpha)/2}) \le C_{2\alpha} \left( \frac{|U|}{|V|} \right)^{\alpha}, $$

where the constant $C_{2\alpha}$ is the one of Lemma F.7. By using Jensen inequality once more, we deduce that

$$ \mathbb{E} \left[ \|W - p\epsilon |V| N^{-1} 1_U \|_2^\alpha \right] \le \left( \mathbb{E} \left[ \|W - p\epsilon |V| N^{-1} 1_U \|_2^{2\alpha} \right] \right)^{1/2}, $$

and the result follows, putting $K_\alpha := C_{2\alpha}^{1/2}$. $\square$

**LEMMA F.9.** Consider a vector $v = (v_1, \dots, v_N) \in \mathbb{R}^N$ supported on $S \subseteq \{1, \dots, N\}$. For any $\xi \in \mathbb{R}$, the following inequality holds:

$$ (95) \qquad \left\| v - \frac{N}{|S|} \bar{v} 1_S \right\|_2 \le \|v - \xi 1_S\|_2, $$

where $\bar{v} = N^{-1} \sum_{i=1}^N v_i = N^{-1} \sum_{i \in S} v_i$.
PROOF. Consider $U \sim \text{Unif}(\{v_i : i \in S\})$, and observe that $\mathbb{IE}[U] = \bar{v}N|S|^{-1}$ and $\text{Var}(U) = |S|^{-1}\|v - N|S|^{-1}\bar{v}1_S\|_2$. Since $\text{Var}(U) = \inf_{\xi \in \mathbb{R}} \mathbb{IE}[(U - \xi)^2]$, it follows that for any $\xi \in \mathbb{R}$,

$$|S|^{-1} \left\| v - \frac{N}{|S|} \bar{v} 1_S \right\|_2 \le \mathbb{IE}[(U - \xi)^2] = |S|^{-1} \|v - \xi 1_S\|_2,$$

implying the result.

<figure><img src="image_16.png" /></figure>

PROOF OF LEMMA F.1. Inequalities (81) and (82) follow from Lemmas F.9 and F.8 respectively. Next, we prove Inequality (80). First, note that $(r_a^N r_b^N)^{-1} \overline{L^{N,ab}} - (bp) = b(|\mathcal{P}_a||\mathcal{P}_b|)^{-1} \sum_{i \in \mathcal{P}_a} \sum_{j \in \mathcal{P}_b} (\theta_{ij} - p)$, so that

$$ (96) \qquad \begin{aligned} \mathbb{IE} \left[ \left| (r_a^N r_b^N)^{-1} \overline{L^{N,ab}} - (bp) \right|^2 \right] &= \text{Var} \left( (|\mathcal{P}_a||\mathcal{P}_b|)^{-1} \sum_{i \in \mathcal{P}_a} \sum_{j \in \mathcal{P}_b} \theta_{ij} \right) \\ &= p(1-p)(|\mathcal{P}_a||\mathcal{P}_b|)^{-1} = N^{-2} p(1-p) (r_a^N r_b^N)^{-1}. \end{aligned} $$

Since $p(1-p) \le 1/4$, we conclude that

$$\mathbb{IE} \left[ \left| (r_a^N r_b^N)^{-1} \overline{L^{N,ab}} - (bp) \right|^2 \right] \le (N^{-2}/4) \max_{c,d} \{(r_c^N r_d^N)^{-1}\}.$$

Similar arguments can be applied to show that

$$\mathbb{IE} \left[ \left| (r_a^N r_b^N)^{-1} \overline{C^{N,ab}} - (ap) \right|^2 \right] \le (N^{-2}/4) \max_{c,d} \{(r_c^N r_d^N)^{-1}\}.$$

Inequality (80) follows then from the above inequalities and from the fact that $x \vee y \le x + y$ for any $x, y \ge 0$.

We now establish Inequality (79). To that end, first observe that

$$\mathbb{IE} \left[ \| A^N \tilde{L}^{N,ab} \|_2^2 \right] = N \mathbb{IE} \left[ \left( \sum_{j \in \mathcal{P}_a} A^N(1,j) \tilde{L}_j^{N,ab} \right)^2 \right].$$

Next, observe that if we denote $Z_j = L_j^{N,ab} - (pb)r_b^N$ for $j \in \mathcal{P}_a$, then we can write (recall that $A^N(i,j) = bN^{-1}\theta_{ij}$ for $1 \le i \le N$ and $j \in \mathcal{P}_b$),

$$\left[ \sum_{j \in \mathcal{P}_a} A^N(1,j) \tilde{L}_j^{N,ab} \right]^2 = N^{-2} \left[ \sum_{j \in \mathcal{P}_a} \theta_{1j} Z_j - \left( \sum_{j \in \mathcal{P}_a} \theta_{1j} \right) \left( |\mathcal{P}_a|^{-1} \sum_{k \in \mathcal{P}_a} Z_k \right) \right]^2,$$

so that by applying Jensen's inequality we obtain that

$$ (97) \qquad \begin{aligned} \mathbb{IE} \left[ \| A^N \tilde{L}^{N,ab} \|_2^2 \right] \le \\ 2N^{-1} \left[ \mathbb{IE} \left[ \left( \sum_{j \in \mathcal{P}_a} \theta_{1j} Z_j \right)^2 \right] + \mathbb{IE} \left[ \left( \sum_{j \in \mathcal{P}_a} \theta_{1j} \right)^2 \left( |\mathcal{P}_a|^{-1} \sum_{k \in \mathcal{P}_a} Z_k \right)^2 \right] \right]. \end{aligned} $$
By observing that $\sum_{j \in \mathcal{P}_a} \theta_{1j} \le |\mathcal{P}_a|$ and using that $(Z_k)_{k \in \mathcal{P}_b}$ are i.i.d. centered random variables, one can check that

$$ \begin{aligned} \mathbb{E} \left[ \left( \sum_{j \in \mathcal{P}_a} \theta_{1j} \right)^2 \left( |\mathcal{P}_a|^{-1} \sum_{k \in \mathcal{P}_a} Z_k \right)^2 \right] &\le |\mathcal{P}_a|^2 \mathrm{Var} \left( |\mathcal{P}_a|^{-1} \sum_{k \in \mathcal{P}_a} Z_k \right) \\ &= |\mathcal{P}_a| \mathrm{Var} (Z_1) = r_a^N r_b^N p(1-p) \le 1/4, \end{aligned} $$

where in the last inequality we have used that $\mathrm{Var}(Z_k) = \mathrm{Var}(L_k^{N,ab}) = |\mathcal{P}_b|p(1-p)/N^2$.

Now, note that if $j \ne 1$, then $\theta_{1j}$ is independent of $Z_j$ so that $\mathbb{E}[\theta_{1j} Z_j] = \mathbb{E}[\theta_{1j}] \mathbb{E}[Z_j] = 0$, because $Z_j$ is centered. Moreover, one can check that $(\theta_{1j} Z_j)_{j \in \mathcal{P}_a : j \ne 1}$ are independent. Hence, by combining these facts with Jensen's inequality, we can deduce that

$$ \begin{aligned} \mathbb{E} \left[ \left( \sum_{j \in \mathcal{P}_a} \theta_{1j} Z_j \right)^2 \right] &\le 2 \left[ \mathbb{E} \left[ \left( \sum_{j \in \mathcal{P}_a : j \ne 1} \theta_{1j} Z_j \right)^2 \right] + \mathbb{E} \left[ (\theta_{11} Z_1)^2 \right] \right] \\ &\le 2 \left[ \mathrm{Var} \left( \sum_{j \in \mathcal{P}_a : j \ne 1} \theta_{1j} Z_j \right) + \mathrm{Var} (Z_1) \right] \\ &\le 2 \left[ (|\mathcal{P}_a| - 1)p \mathrm{Var}(Z_1) + \mathrm{Var}(Z_1) \right] \le 2 r_a^N r_b^N p(1-p) \le 1/2. \end{aligned} $$

Combining the last three inequalities we show that Inequality (79) holds.

It remains to show Inequality (78). In the case $a_1 \ne a$, the two vectors $\tilde{L}^{N,ab}$ and $\tilde{L}^{N,ab_1}$ have disjoint supports so that $\langle \tilde{L}^{N,ab}, \tilde{L}^{N,ab_1} \rangle = 0$ a.s. and (79) trivially holds.

From now on, consider that $a_1 = a$. Let us denote $Y_{a,bb_1}^N = \langle \tilde{L}^{N,ab}, \tilde{L}^{N,ab_1} \rangle$ and prove that there exists a universal constant $K$ such that

$$ (98) \quad \mathbb{E} \left[ \left| Y_{a,bb_1}^N - \delta_{bb_1} r_a r_b p(1-p) \right|^2 \right] \le K \left[ \mathbb{E} \left[ \left| Y_{a,bb_1}^N - \mathbb{E} [Y_{a,bb_1}^N] \right|^2 \right] + N^{-2} \right]. $$

First, consider the case $b_1 \ne b$. Remind that $L^{N,ab}$ only depends on $\{\theta_{ij}, i \in \mathcal{P}_a, j \in \mathcal{P}_b\}$, and remark that $\mathbb{E}[\tilde{L}_i^{N,ab}] = 0$ for all $i = 1, \dots, N$. In particular, the vectors $L^{N,ab}$ and $L^{N,ab_1}$ are independent and

$$ \mathbb{E} [Y_{a,bb_1}^N] = \langle \mathbb{E} [\tilde{L}^{N,ab}], \mathbb{E} [\tilde{L}^{N,ab_1}] \rangle = 0, $$

so that

$$ \mathbb{E} \left[ \left| Y_{a,bb_1}^N - \delta_{bb_1} r_a r_b p(1-p) \right|^2 \right] = \mathbb{E} \left[ \left| Y_{a,bb_1}^N - \mathbb{E} [Y_{a,bb_1}^N] \right|^2 \right]. $$

Now, consider the case $b_1 = b$. In that case, $Y_{a,bb_1}^N = \|\tilde{L}^{N,ab}\|_2^2$, and observe that

$$ \mathbb{E} \left[ \| \tilde{L}^{N,ab} \|_2^2 \right] = |\mathcal{P}_a| \mathbb{E} \left[ \left( L_i^{N,ab} - \frac{1}{|\mathcal{P}_a|} \sum_{j \in \mathcal{P}_a} L_j^{N,ab} \right)^2 \right], $$

for any fixed $i \in \mathcal{P}_a$. In the rest of the proof, $i$ denotes an arbitrary index in $\mathcal{P}_a$. Next, observe that if we denote $Z_j = L_j^{N,ab} - (pb)r_b^N$ for $j \in \mathcal{P}_a$, then we can write

$$ \mathbb{E} \left[ \| \tilde{L}^{N,ab} \|_2^2 \right] = |\mathcal{P}_a| \mathbb{E} \left[ \left( Z_i(1 - |\mathcal{P}_a|^{-1}) - \frac{1}{|\mathcal{P}_a|} \sum_{j \in \mathcal{P}_a : j \ne i} Z_j \right)^2 \right]. $$