# Theory and inference for multivariate autoregressive binary models with an application to absence-presence data in ecology

Guillaume Franchi, Lionel Truquet

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

## To cite this version:

Guillaume Franchi, Lionel Truquet. Theory and inference for multivariate autoregressive binary models with an application to absence-presence data in ecology. 2025. hal-04996072

HAL Id: hal-04996072

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

Preprint submitted on 18 Mar 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.

Public Domain
# Theory and inference for multivariate autoregressive binary models with an application to absence-presence data in ecology

Guillaume Franchi guillaume.franchi@ensai.fr

Univ Rennes, Ensai, CNRS, CREST – UMR 9194, F-35000 Rennes, France

Lionel Truquet lionel.truquet@ensai.fr

Univ Rennes, Ensai, CNRS, CREST – UMR 9194, F-35000 Rennes, France

## Abstract

We introduce a general class of autoregressive models for studying the dynamic of multivariate binary time series with stationary exogenous covariates. Using a high-level set of assumptions, we show that existence of a stationary path for such models is almost automatic and does not require parameter restrictions when the noise term is not compactly supported. We then study in details statistical inference in a dynamic version of a multivariate probit type model, as a particular case of our general construction. To avoid a complex likelihood optimization, we combine pseudo-likelihood and pairwise likelihood methods for which asymptotic results are obtained for a single path analysis and also for panel data, using ergodic theorems for multi-indexed partial sums. The latter scenario is particularly important for analyzing absence-presence of species in Ecology, a field where data are often collected from surveys at various locations. Our results also give a theoretical background for such models which are often used by the practitioners but without a probabilistic framework.

**Keywords** — Autoregressive Probit Time Series; Pairwise Likelihood Estimation; Panel Data; M-Estimators

## 1 Introduction

As a benchmark model for studying several correlated binary responses, multivariate logit or probit regression models have been extensively used in various domains. For instance, they have been applied to labor supply in economics (Chib & Greenberg, 1998), to the modeling of multiple claims in insurance (Young et al., 2009), for evaluating customer satisfaction (Kenne Pagui & Canale, 2016), for studying the response of biological organisms exposed to stimulus (Ashford & Sowden, 1970) or for studying absence-presence of several species in an ecosystem (Wilkinson et al., 2019), (Poggiato et al., 2021). Sebastian-Gonzalez et al. (2010) also applied a multivariate logistic model for analyzing such absence-presence time series in several ecosystems and Guanche et al. (2014) studied atmospheric circulation with similar models.

Though a few contributions considered dynamic multivariate probit models for time series analysis of binary vectors (see for instance Candelon et al. (2013), Chaubert et al. (2008), Manner et al. (2016)), theoretical contributions to statistical inference in such models are rather scarce. A theory for the univariate dynamic probit model can be found in De Jong and Woutersen (2011) while Tuzcuoglu (2023) investigated the AR-probit model introduced by Müller and Czado (2005) in the case of panel data.
Let us also mention another type of multiple binary autoregressive time series models based on random binary coefficients, an approach developed recently by Jentsch and Reichmann (2022) who introduced a rigorous framework for defining stationary solutions and estimating model parameters. This class of model seems to be quite flexible but a possible inclusion of exogenous covariates in the corresponding model equations is unclear because the link function cannot map continuous data to multiple binary vectors.

In this paper, we consider the following model. For $i \in \{1, \dots, k\}$, let $(Y_{i,t})_{t \in \mathbb{Z}}$ be a univariate binary time series. Following our initial motivation for applying these models to Ecology, we will often refer to a specific application in this context, where $Y_{i,t}$ simply measures the absence-presence of a specific species $i$ at time $t$ in a given area. This analogy will be often useful to interpret model parameters, though our generic results can be used for analyzing arbitrary multiple binary time series. Setting $Y_t = (Y_{1,t}, \dots, Y_{k,t}) \in \{0, 1\}^k$, we assume that the dynamic of the multivariate time series $(Y_t)_{t \in \mathbb{Z}}$ is given by

$$Y_{i,t} = 1_{(0, +\infty)} (g_i(Y_{t-1}, \dots, Y_{t-p}, \zeta_t)), \quad 1 \le i \le k, \quad (1)$$

where $g = (g_1, \dots, g_k)' : \{0, 1\}^{kp} \times \mathbb{R}^n \to \mathbb{R}^k$ is a measurable mapping and $(\zeta_t)_{t \in \mathbb{Z}}$ is a stationary stochastic process taking values in $\mathbb{R}^n$.

Hence the multivariate binary vector $Y_t$ depends on its values at previous times and can be interpreted as the discretization of a multivariate continuous response $Y_t^* = g(Y_{t-1}, \dots, Y_{t-p}, \zeta_t)$. Allowing dependence between the variables $(\zeta_t)_{t \in \mathbb{Z}}$ is interesting for including exogenous covariates, i.e. $\zeta_t = (X_{t-1}, \varepsilon_t)$, where $X_{t-1}$ is a vector of covariates determined at previous time (and then observable before time $t$) and $\varepsilon_t$ is an unobserved noise component (that in some cases will be assumed to be independent from $X_{t-1}$). When $g$ is linear, i.e.

$$g(Y_{t-1}, \dots, Y_{t-p}, \zeta_t) = \sum_{j=1}^{p} A_j Y_{t-j} + B X_{t-1} + \varepsilon_t, \quad (2)$$

where $A_1, \dots, A_p$ are square matrices of size $k \times k$, $B$ is a matrix of size $k \times d$, $\varepsilon_t$ is independent of $(Y_{t-1}, \dots, Y_{t-p}, X_{t-1})$ and follows a Gaussian distribution with mean 0 and covariance matrix $R$, we obtain a dynamic version of the multivariate probit regression model.

More generally, model (1) is a natural time series version of the multivariate discrete choice models popular in econometrics (see for instance Greene (2009)). For $k = 1$, it extends the dynamic binary choice model of De Jong and Woutersen (2011). When $F$ is the CDF of the logistic distribution, we get the multivariate logistic model used for instance by Sebastian-Gonzalez et al. (2010) for studying absence-presence data. Note that our version has some natural interpretations in term of dynamics for a multivariate binary response. For instance, when $p = 1$ for simplicity, $A_1(i, j)$ is positive if and only if the presence of a species $j$ at time $t-1$ influences positively the presence of the species $i$ at time $t$. Furthermore, the sign of $R_{i,j}$ indicates if positive or negative dependence occurs at time $t$ (positivity or negativity of the conditional covariances between $Y_{i,t}$ and $Y_{j,t}$ meaning that the values 1 are generally produced simultaneously or not).

This paper provides a threefold contribution.

First, we give simple and very general conditions guarantying the existence of a strictly stationary solution for (1). This point, which is not straightforward because of the non-linearity of the link function, is always an important step in time series analysis since it is very useful to derive asymptotic properties of statistical inference procedures. As we will see, stationarity is almost automatic for model (2), without any restriction on the coefficients. Note that from the presence of exogenous covariates, the model is not Markovian in general, leading to an additional difficulty for studying stationarity properties.

Secondly, we investigate statistical inference for the autoregressive parameters $A_1, \dots, A_p, B$ and the correlation matrix $R$ in model (2), which is quite challenging. Likelihood estimators are generally intractable due to the presence of multiple integrals involving the Gaussian density in the expression of the likelihood function. There exists many approaches developed for inference in multivariate probit
models. Some of them are based on a Bayesian framework (Chib & Greenberg, 1998), (Talhouk et al., 2012), or on importance sampling such as the GHK algorithm (Hajivassiliou & Ruud, 1994), or on composite likelihood methods (Zhao & Joe, 2005), (Ting et al., 2022). In the time series context, we will combine univariate and pairwise composite likelihood methods to estimate all the parameters of the model. The advantage of the proposed method is its fast implementation with respect to the methods based on the true likelihood function. We are not aware of a theoretical analysis of some inference procedures for multivariate logistic/probit models such as the double-step approach mentioned above.

Finally, the case of data collected at different sites is even more challenging due to a double asymptotic when both the number of sites $n$ and the horizon time $T$ go to infinity. The benefit of this work is to provide very simple assumptions for getting asymptotic results in this context. In particular, we use specific versions of the law of large numbers for partial sums with a double summation index. All the results are formulated using realistic assumptions on the covariates which can contain individual effects, time effects or both. Let us also mention that our results, which are based on a multi-parameter ergodic theorem not very well known in the time series literature, can be also interesting for studying more general panel type models.

The paper is organized as follows. In Section 2, we address the problem of stationarity for these dynamics. Inference methods for the multivariate autoregressive probit model are introduced and discussed for a single path in Section 3, while Section 4 is devoted to the problem of multiple paths with panel data. The results are illustrated in Section 5 with an application to the dynamics of fisheries data. The proofs of all our results are given in a supplementary material available online.

# 2 Existence of stationary solutions

## 2.1 A general result for multiple binary autoregressive model

In this section, we first state a general result ensuring the existence of a stationary solution for model (1). The assumptions used below will be much more general than what is needed for studying the multivariate linear model (2). In what follows, we denote by $\mathcal{F}_t$ the sigma-field generated by the random variables $\zeta_s$, $s \le t$:

$$ \mathcal{F}_t = \sigma(\zeta_s | s \le t). $$

We consider the following set of assumptions.

**A1** The process $(\zeta_t)_{t \in \mathbb{Z}}$, taking values in $\mathbb{R}^n$, is strictly stationary.

**A2** For $i = 1, \dots, k$, let us introduce the two events

$$ C_i^0 = \left\{ \max_{1 \le t \le p} \max_{y_1, \dots, y_p \in \{0,1\}^k} g_i(y_1, \dots, y_p, \zeta_t) \le 0 \right\}, $$

and

$$ C_i^1 = \left\{ \min_{1 \le t \le p} \min_{y_1, \dots, y_p \in \{0,1\}^k} g_i(y_1, \dots, y_p, \zeta_t) > 0 \right\}. $$

We assume that there exist $j_1, \dots, j_k$ in $\{0, 1\}$ such that $\mathbb{P}(\cap_{i=1}^k C_i^{j_i}) > 0$.

**A3** There exist $j_1, \dots, j_k$ in $\{0, 1\}$ such that $\mathbb{P}(\cap_{i=1}^k C_i^{j_i} | \mathcal{F}_0) > 0$ almost surely.

Let us first comment these assumptions. Assumption **A1** is the minimal assumption needed for working with stationary time series. The technical assumption **A2** has an important interpretation. Assume that the stochastic recursions in (1) are initialized with some values $Y_t = y_t$ for $-p+1 \le t \le 0$. If the event $\cap_{i=1}^k C_i^{j_i}$ occurs, whatever the initial values $y_{-p+1}, \dots, y_0$, we will obtain the value $Y_{i,t} = j_i$.
for $1 \le t \le p$ and $1 \le i \le k$. As a consequence, the values of $Y_t$ for $t \ge p+1$ will no more depend on the initial values. This kind of condition is analog to the standard irreducibility and aperiodicity conditions imposed to finite-state Markov chains though the model studied here does not have the Markov property.

From these two assumptions, it is possible to prove a stochastic convergence property for the iterations of the recursions (1). The proof of the following result is based on a general coalescence argument for iterated random maps acting on finite sets, an idea already used by Z. M. Debaly and Truquet (2021) for the binary choice model and extended by Truquet (2023) to more general model, in the case of ergodic inputs. This argument has some flavors with a similar one used by Propp and Wilson (1996) to simulate perfectly a realization of the invariant probability of a finite-state Markov chain. However, our aim here is also to discuss the existence of stationary but not necessarily ergodic paths (as it will be required for studying panel data with individual effects). We have to adapt the techniques used in these references, replacing Assumption A2 by a more stringent condition given in Assumption A3.

**Theorem 2.1.**

1) Suppose that assumptions **A1-A2** hold true and the process $(\zeta_t)_{t \in \mathbb{Z}}$ is ergodic. Then there exists a unique strictly stationary solution $(Y_t)_{t \in \mathbb{Z}}$ for the recursions (1). Moreover, the unique solution has a Bernoulli shift representation, i.e. there exists a measurable mapping $H : (\mathbb{R}^n)^{\mathbb{N}} \to \{0, 1\}$ such that

$$Y_t = H(\zeta_t, \zeta_{t-1}, \dots), \quad t \in \mathbb{Z}. \qquad (3)$$

The process $((Y_t, \zeta_t))_{t \in \mathbb{Z}}$ is then ergodic.

2) Suppose now that assumptions **A1-A3** hold true with $(\zeta_t)_{t \in \mathbb{Z}}$ not necessarily ergodic. There still exists a unique strictly stationary solution $(Y_t)_{t \in \mathbb{Z}}$ for the recursions (1) and the Bernoulli shift representation is still valid.

**Notes.**

1. When $\zeta_t = (X_{t-1}, \varepsilon_t)$ and $X_t$ is a vector of covariates observed at time $t$ and $\varepsilon_t$ is an unobserved noise component, we point out that the non-ergodic case is important in some applications. Indeed, in panel data, some of the components of the vector $X_t$ may be not time-varying. For instance, $X_t = (Z, \Gamma_t)$ with $Z$ being a random vector independent of a stationary process $((\Gamma_t, \varepsilon_t))_{t \in \mathbb{Z}}$. In this case, the process $(\zeta_t)_{t \in \mathbb{Z}}$ is still stationary but obviously not ergodic. Let us also mention that stationarity of the process $(\zeta_t)_t$ is equivalent to stationarity for the process $(X_t, \varepsilon_t)_t$ and a similar remark is valid for ergodicity.

2. In the ergodic case, Theorem 2.1, point 1) is useful for applying Birkhoff's ergodic theorem to some partial sums of the form $\frac{1}{T} \sum_{t=1}^{T} f(Y_t, X_t, Y_{t-1}, X_{t-1}, \dots)$ where $f$ is a measurable function such that $\mathbb{E}[f(Y_0, X_0, Y_{-1}, X_{-1}, \dots)] < +\infty$. This ergodic theorem is generally sufficient to derive consistency and asymptotic normality of conditional likelihood or pseudo-likelihood estimators in smooth parametric time series models.

## 2.2 Specific results for models with covariates and an additive noise component

In this part, we assume that $\zeta_t = (X_{t-1}, \varepsilon_t)$ takes values in $\mathbb{R}^d \times \mathbb{R}^k$ (then $n = d+k$) and we consider the more specific dynamic

$$g(Y_{t-1}, \dots, Y_{t-p}, \zeta_t) = h(Y_{t-1}, \dots, Y_{t-p}, X_{t-1}) + \varepsilon_t, \qquad (4)$$

where $h: (\{0, 1\}^k)^p \times \mathbb{R}^d \to \mathbb{R}^k$ is a measurable mapping. Of course, one could also assume a dependence with respect to $X_{t-2}, \dots, X_{t-p}$ but in this case it is possible to replace $X_{t-1}$ in (4) by the vector $Z_{t-1} := (X_{t-1}, \dots, X_{t-p})$ so that model (4) contains implicitly such an extension.
The following condition is sufficient for getting **A3** and then **A2**.

**A4** Almost surely, the conditional distribution of $\varepsilon_t$ given $\sigma(\varepsilon_s, X_s | s \le t - 1)$ has a full support $\mathbb{R}^k$.

Assumption **A4** is sufficient but not necessary for checking **A3**, other more sophisticated conditions on the support of the previous conditional distribution are possible, but **A4** is simpler to formulate. From Assumption **A4**, we get a simple and meaningful condition for recurrence of the states in the dynamic. Indeed, whatever the previous values of the process and the covariates are, the process is allowed to move to an arbitrary element of $\{0, 1\}^k$ (i.e. we have positive conditional probabilities).

Note that Assumption **A4** does not require a specific form for the mapping $h$. Moreover, it is valid when $\varepsilon_t$ is independent from $\sigma(\varepsilon_s, X_s | s \le t - 1)$ (i.e. the covariates are sequentially exogenous) and its probability distribution has a full support (as for the multivariate Gaussian for instance). But assumption **A4** also holds true with non i.i.d. noises. For instance, it is satisfied for infinite moving averages $\varepsilon_t = \eta_t + \sum_{j=1}^{\infty} A_j \eta_{t-j}$ where $(A_j)_{j \ge 0}$ is a summable sequence of square matrices and $((X_t, \eta_t))_{t \in \mathbb{Z}}$ is a stationary process such that $\eta_t$ has a probability distribution with full support and is independent from $\sigma((\eta_s, X_s) | s \le t - 1)$. Let us recall that from the definition of the support of a probability measure, i.e. the intersection of all the closed sets with probability 1, the full support assumption entails that $\mathbb{P}(\eta_t \in O) > 0$ for any non-empty open subset $O$.

**Corollary 2.2.** Suppose that Assumptions **A1-A4** hold true.

Then the conclusions of Theorem 2.1 are valid.

**Note.** When there is no covariates (i.e. $h$ does not depend on $X_{t-1}$) and the $\varepsilon_t$'s are i.i.d., the process is simply an irreducible $p$-order Markov chain and we already know that there exists a unique strictly stationary solution to (1) and (2). Although the process is not necessarily Markovian here due to the presence of the covariates, one can observe that the stationarity conditions are similar and do not impose any specific condition on $h$. Things become more complicated if the noise $\varepsilon_t$ has a bounded support even when no covariates are included in the dynamics. For instance, suppose that $p = k = 1$ with $\varepsilon_t$ supported on $[0, 1]$. When $h(y) = 2y - 1$, the chain is not irreducible and points 0 and 1 are absorbing states. On the other hand, when $h(y) = -y$, the chain is irreducible but periodic as it alternates between state 0 and state 1. One can then see on this simple example that the behavior of the chain strongly depends on the parameters of the model.

# 3 Theory and inference for a multivariate probit autoregressive model

## 3.1 Framework

In this section, we consider a stochastic process $(Y_t)_{t \in \mathbb{Z}}$ satisfying equations (1) and (2) with multivariate Gaussian errors i.e.

$$ \forall t \in \mathbb{Z}, \forall i \in \{1, \dots, k\}, Y_{i,t} = \mathbb{1}_{(0, +\infty)} (\lambda_{i,t} + \varepsilon_{i,t}), \quad (5) $$

where

$$ \lambda_t = C + \sum_{l=1}^{p} A_l Y_{t-l} + B X_{t-1}. \quad (6) $$

The sequences $(X_t)_{t \in \mathbb{Z}}$ and $(\varepsilon_t)_{t \in \mathbb{Z}}$ of random vectors are valued in $\mathbb{R}^d$ and $\mathbb{R}^k$ respectively, $A_1, \dots, A_p$ are matrices of size $k \times k$, $B$ is a matrix of size $k \times d$ and $C$ is a vector of $\mathbb{R}^k$. The process $(X_t)_{t \in \mathbb{Z}}$ is actually a process of covariates, while $(\varepsilon_t)_{t \in \mathbb{Z}}$ corresponds to an unobserved noise component.
Finally, we now denote by $\mathcal{F}_t$ the sigma-field generated by $X_s$ and $\varepsilon_s$ for $s \le t$, i.e.

$$ \mathcal{F}_t = \sigma (X_s, \varepsilon_s | s \le t). $$

The following assumptions will be used.

B1 The process $(\zeta_t)_{t \in \mathbb{Z}}$ is strictly stationary (resp. strictly stationary and ergodic).

B2 For any $t \in \mathbb{Z}$, the random vector $\varepsilon_t$ is independent from $\mathcal{F}_{t-1}$.

B3 The random vectors $\varepsilon_t$ are normally distributed

$$ \varepsilon_t \sim \mathcal{N}(0, R) $$

where $R$ is a correlation matrix.

The following result is a straightforward consequence of Corollary 2.2.

**Corollary 3.1.** Suppose that Assumptions B1-B2-B3 hold true. Then there exists a unique stationary (resp. stationary and ergodic) solution $(Y_t))_{t \in \mathbb{Z}}$ to the recursions (5) and (6).

Notes.

1. Assumptions **B1** and **B2** are verified in the particular case where there exists a measurable application $H$ such that $X_t = H(\eta_t, \eta_{t-1}, \dots)$, with the variables $(\eta_t, \varepsilon_t)_{t \in \mathbb{Z}}$ i.i.d., but $\eta_t$ not necessarily independent from $\varepsilon_t$. For example, $(X_t)_{t \in \mathbb{Z}}$ can be a MA($\infty$) process. This assumption is sufficiently general and does not require a specific dynamic for the covariates such as a Markov property for instance.

2. When $A_1 = \dots = A_p = 0$, the previous model reduces to the standard multivariate probit regression model studied for instance in Chib and Greenberg (1998). We also recover a specific case of the dynamic model studied in Candelon et al. (2013) but without theory. Note that our version is particularly convenient for model interpretation. For instance, in the case of absence/presence data, a positive (resp. negative) sign for $A_j(i, \ell)$ indicates if the presence of an element $\ell$ at time $t-j$ promotes or not the presence of the element $i$ at time $t$. An interpretation of the elements of the correlation matrix $R$ can be also given in term of conditional dependencies. For instance, conditionally on $\mathcal{F}_{t-1}$, the conditional covariance between $Y_{i,t}$ and $Y_{\ell,t}$ have the same sign as $R(i, \ell) \in [-1, 1]$. Indeed if $\Phi$ denotes the cdf of the standard Gaussian distribution and if $\Phi_{R(i,\ell)}$ denotes the cdf of a Gaussian vector with standard Gaussian marginals and correlation coefficient $R(i, \ell)$, we have

$$ \mathrm{Cov}(Y_{i,t}, Y_{\ell,t} | \mathcal{F}_{t-1}) = \Phi_{R(i,\ell)}(\lambda_{i,t}, \lambda_{\ell,t}) - \Phi(\lambda_{i,t}) \Phi(\lambda_{\ell,t}). $$

From Slepian's lemma, this quantity increases with $R(i, \ell)$ and equals to 0 when $R(i, \ell) = 0$.

3. Another very nice property of this model is its invariance if we are only interested in predicting or explaining a specific subset of the $k$ binary random variables. Indeed, if equation (6) is valid for the $k$ coordinates, then for any subset $I$ of $\{1, \dots, k\}$, the time series $((Y_{i,t})_{i \in I})_{t \in \mathbb{Z}}$ follows the same kind of dynamic by considering the $Y'_{i,t-j}$s for $i \notin I$ as exogenous regressors as we do for $X_{t-1}$. Including exogenous regressors satisfying **B2**, without assuming complete independence with the entire noise sequence is then very important for such models.

4. Considering a correlation matrix $R$ instead of a covariance matrix for the noise term $\varepsilon_t$ is mainly motivated by identification issues. One can always multiply $\lambda_t$ and $\varepsilon_t$ by a diagonal matrix with positive entries and this will not change the sign of the coordinates. Up to a change of parameters, one can always assume the covariance matrix of $\varepsilon_t$ to be a correlation matrix.
## 3.2 Statistical inference from a single path

In this part, we assume that the random vectors $Y_1, \dots, Y_T$ and $X_1, \dots, X_{T-1}$ are observed for some positive integer $T$. Our aim is to estimate regression parameters $C, A_1, \dots, A_p, B$ but also the correlation matrix $R$ of the noise component $\varepsilon_t$. Likelihood inference in multivariate probit model is particularly tricky due to the presence of multiple integrals of the Gaussian density over some rectangles with boundaries depending on the unknown parameters. Bayesian inference combined with MCMC methods have been developed in Chib and Greenberg (1998) and there also exist importance sampling strategies such as in Liesenfeld and Richard (2010) or Masarotto and Varin (2012). In the dynamic setting, the conditional likelihood is based on the conditional distribution given below. We first introduce some notations. We set

$$I_1 = (0, +\infty) \quad \text{and} \quad I_0 = (-\infty, 0].$$

The vector of all unknown parameters is given by

$$\theta = (\gamma, r')$$

where $\gamma$ is a column vector containing the column vectors of the matrices $C, A_1, \dots, A_p, B$ and

$$r = (R(1, 2), \dots, R(1, k), R(2, 3), \dots, R(2, k), \dots, R(k-1, k)).$$

The vector $\lambda_t$ is actually a function of $\gamma$, and without ambiguity, we will rather denote $\lambda_t$ by $\lambda_t(\gamma)$ and $\lambda_{i,t}$ by $\lambda_{i,t}(\gamma)$ for $1 \le i \le k$.

In what follows, we focus on two estimation procedures that can be used for the inference of the true parameter $\theta_0 = (\gamma_0, r'_0)$ satisfying (5) and (6). For $t \in \mathbb{Z}$, the distribution of $Y_t$ conditionally to $\mathcal{F}_{t-1}$ is given by

$$\mathcal{L}(Y_t | \mathcal{F}_{t-1}) = \sum_{s \in \{0,1\}^k} \left( \int_{\mathbb{R}^k} \prod_{i=1}^k 1_{I_{s(i)}} (\lambda_{i,t}(\gamma_0) + x_i) \varphi_{R_0}(x) dx \right) \delta_s,$$

where $\varphi_R$ denotes the Gaussian density on $\mathbb{R}^k$ with mean 0 and covariance matrix $R$. When $k=1$, we simply denote $\varphi_R$ by $\varphi$.

Therefore we define the pseudo conditional log-likelihood function by

$$\mathcal{PL}(\theta) = \sum_{t=p+1}^{T} \log \left( \int_{\mathbb{R}^k} \prod_{i=1}^{k} 1_{I_{Y_{i,t}}} (\lambda_{i,t}(\gamma) + x_i) \varphi_R(x) dx \right), \quad (7)$$

and we set $\hat{\theta}_{\mathcal{PL}} = \underset{\theta \in \Theta}{\operatorname{argmax}} \mathcal{PL}(\theta)$.

Note that (7) is based on the density of $(Y_t)_{p+1 \le t \le T}$ conditionally on $Y_1, \dots, Y_p$ and $X_1, \dots, X_T$. However, such a conditional density coincides with (7) only if the $X_t$'s are independent from the $\varepsilon_t$'s, which is not necessarily the case in our framework, hence the terminology pseudo conditional log-likelihood. One main drawback of (7) is the computation of many iterated integrals with some random bounds depending on the autoregressive parameters. Hence, numerical evaluations of the conditional likelihood function is complicated.

We will thus study a two-step method by first estimating the autoregressive parameters $A_1, \dots, A_p, B$ and then the correlation matrix $R$. We first define for a positive integer $T \ge p+1$

$$\begin{aligned} \overline{\mathcal{PL}}(\gamma) &= \sum_{t=p+1}^{T} \left\{ \sum_{i=1}^{k} Y_{i,t} \log [1 - \Phi(-\lambda_{i,t}(\gamma))] + (1 - Y_{i,t}) \log [\Phi(-\lambda_{i,t}(\gamma))] \right\} \\ &= \sum_{t=p+1}^{T} \left\{ \sum_{i=1}^{k} Y_{i,t} \log [\Phi(\lambda_{i,t}(\gamma))] + (1 - Y_{i,t}) \log [\Phi(-\lambda_{i,t}(\gamma))] \right\}, \end{aligned} \quad (8)$$

where $\Phi$ denotes the CDF of the standard Gaussian distribution in $\mathbb{R}$.
We also set $\hat{\gamma} = \underset{\gamma \in \Gamma}{\operatorname{argmax}} \overline{\mathcal{PL}}(\gamma)$.

Note that (8) is simply a pseudo conditional log-likelihood function. It is obtained from (7) under an additional independence assumption on the coordinates of $\varepsilon_t$ for any $t \in \mathbb{Z}$.

Note also that the maximization of (8) can be obtained equation by equation. Indeed, we have for $1 \le i \le k$ and $t \in \mathbb{Z}$

$$ \lambda_{i,t}(\gamma) = \sum_{l=1}^{p} \sum_{j=1}^{k} A_l(i,j)Y_{j,t-l} + \sum_{j=1}^{d} B(i,j)X_{j,t-1} $$

and only the lines number $i$ of the matrices are involved in the previous expression. Hence maximizing (8) is equivalent to maximize $k$ functions with different parameters.

Once we have an estimator $\hat{\gamma}$ obtained by optimizing (8), one can plug it in the expression of the conditional likelihood function (7) given in the previous point, and maximize it with respect to $R$. However, this still requires numerical evaluations of multiple integrals.

We prefer to adopt a pairwise conditional likelihood approach similar to the one presented in de Leon (2005) (see Varin et al. (2011) for a review about composite likelihood approaches). For $i < j$, we define

$$ R_{i,j} = \begin{pmatrix} 1 & r_{i,j} \\ r_{i,j} & 1 \end{pmatrix} \quad \text{where } r_{i,j} = R(i,j), $$

and

$$ \begin{aligned} \hat{r}_{i,j} &= \underset{r_{i,j} \in \mathcal{R}_{i,j}}{\operatorname{argmax}} \sum_{t=p+1}^{T} \log \left( \int_{I_{Y_{i,t}-\lambda_{i,t}(\hat{\gamma})}} \int_{I_{Y_{j,t}-\lambda_{j,t}(\hat{\gamma})}} \varphi_{R_{i,j}}(x_i, x_j) dx_i dx_j \right) \\ &= \underset{r_{i,j} \in \mathcal{R}_{i,j}}{\operatorname{argmax}} \sum_{t=p+1}^{T} \log \left\{ \int_{I_{Y_{i,t}-\lambda_{i,t}(\hat{\gamma})}} \Phi \left( (2Y_{j,t} - 1) \frac{\lambda_{j,t}(\hat{\gamma}) + r_{i,j}x_i}{\sqrt{1 - r_{i,j}^2}} \right) \varphi(x_i) dx_i \right\}. \end{aligned} \quad (9) $$

Note that for any $t \in \mathbb{Z}$, and $s_i, s_j \in \{0, 1\}$,

$$ \int_{I_{s_i - \lambda_{i,t}(\gamma_0)}} \Phi \left( (2s_j - 1) \frac{\lambda_{j,t}(\gamma_0) + r_{i,j}x_i}{\sqrt{1 - r_{i,j}^2}} \right) \varphi(x_i) dx_i = \mathbb{P} (Y_{i,t} = s_i, Y_{j,t} = s_j | \mathcal{F}_{t-1}), $$

which explains the terminology pairwise conditional likelihood. However, the optimization (9) does not correspond to the optimization of a likelihood function. Here, we simply sum the log-densities of a given pair $(Y_{i,t}, Y_{j,t})$ along time $t$, and maximize it.

Finally, we denote $\hat{r} = (\hat{r}_{1,2}, \dots, \hat{r}_{1,k}, \hat{r}_{2,3}, \dots, \hat{r}_{2,k}, \dots \hat{r}_{k-1,k})$, and for simplicity of further notations, we set for $t \in \mathbb{Z}$ and $\theta = (\gamma, r') \in \Theta$

$$ \ell_t(\theta) = \sum_{1 \le i < j \le k} \log \left( \int_{I_{s_i - \lambda_{i,t}(\gamma)}} \Phi \left( (2Y_{j,t} - 1) \frac{\lambda_{j,t}(\gamma) + r_{i,j}x_i}{\sqrt{1 - r_{i,j}^2}} \right) \varphi(x_i) dx_i \right). \quad (10) $$

For deriving asymptotic properties of these M-estimators, the following additional assumptions will be needed.

**B4** All matrices $A_1, \dots, A_p$ belong to a compact set $\mathcal{A} \subset \mathbb{R}^{k \times k}$, $B$ belongs to a compact set $\mathcal{B} \subset \mathbb{R}^{k \times d}$ and $C$ belongs to compact set of $\mathbb{R}^k$. In addition, $R$ belongs to a compact set $\mathcal{R}$ included in the set of positive definite symmetric matrices in $\mathbb{R}^{k \times k}$.

**B5** For all $t \in \mathbb{Z}$ and $v \in \mathbb{R}^d \setminus \{0\}$, we have for any $\sigma(\varepsilon_t) \lor \mathcal{F}_{t-1}$-measurable application $G$

$$ \mathbb{P} (v' \cdot X_t = G) < 1 $$
Notes.

1. To construct a compact subset of the set of correlation matrices, one can use Cholesky decomposition $R = LL^T$ with $L$ lower triangular and $L(i, i) = \sqrt{1 - \sum_{\ell=1}^{i-1} L(i, \ell)^2}$. One can then assume that $|L(i, \ell)| \le \bar{L}(i, \ell)$ where the $\bar{L}(i, \ell)$'s satisfy $\sum_{\ell=1}^{i-1} \bar{L}(i, \ell)^2 < 1$ for $i = 1, \dots, k$ and $1 \le \ell < i$. The number of free parameters is $(k^2 - k)/2$. Cholesky decomposition will be used in our numerical experiments.

2. Assumption **B5** guarantees that any non-trivial linear combination of the stochastic covariates at time $t$ cannot be almost surely explained by the past values and the sigma-field $\sigma(\varepsilon_t)$. Note that the error term is not necessarily independent of $X_t$, i.e. contemporaneous dependence is allowed between both quantities. By stationarity, it is sufficient to verify assumption **B5** only for $t = 0$. Under the more restrictive assumption of complete independence between the two processes $(X_t)_{t \in \mathbb{Z}}$ and $(\varepsilon_t)_{t \in \mathbb{Z}}$, the identifiability assumption **B5** can be replaced by the following one which is the standard assumption used for regression models in the i.i.d. setting. Details are given in the proof of Theorem 3.2.

**B5'** For any $t \in \mathbb{Z}$, the distribution of $X_t$ is not supported by an affine hyperplane.

We now derive asymptotic results for the estimators defined above. In what follows, for an arbitrary integer $n \ge 2$, we denote by $\nabla f$ the gradient of a real-valued mapping $f$ defined on a subset of $\mathbb{R}^n$. If $n = 1$, we denote by $\dot{f}$ the derivative of the mapping $f$. If $f$ is defined on $\Theta$, we will denote by $\nabla_1 f(\gamma, r)$ and $\nabla_2 f(\gamma, r)$ the two gradients of the functions $\gamma \mapsto f(\gamma, r)$ and $r \mapsto f(\gamma, r)$ evaluated at point $\theta = (\gamma, r)$.

Similarly, we denote by $\nabla^2 f$ the hessian matrix of $f$, and if $n = 1$, we simply denote by $\ddot{f}$ the second derivative of the mapping $f$. If $f$ is defined on $\Theta$, setting $d_\Gamma$ and $d_\mathcal{R}$ the dimensions of $\Gamma$ and $\mathcal{R}$ respectively, we will denote

$$ \nabla_1^2 f = \left( \frac{\partial^2 f}{\partial \gamma_i \partial \gamma_j} \right)_{1 \le i,j \le d_\Gamma}, \quad \nabla_{1,2}^2 f = \left( \frac{\partial^2 f}{\partial \gamma_j \partial r_i} \right)_{\substack{1 \le i \le d_\mathcal{R} \\ 1 \le j \le d_\Gamma}}, \quad \text{and} \quad \nabla_2^2 f = \left( \frac{\partial^2 f}{\partial r_i \partial r_j} \right)_{1 \le i,j \le d_\mathcal{R}}. $$

Finally, we denote by $\| \cdot \|_2$ the Euclidean norm.

The following theorem establishes the consistency and asymptotic normality of the estimators $\hat{\theta}_{\mathcal{P}\mathcal{L}}$, $\hat{\gamma}$ and $\hat{r}$ defined previously.

**Theorem 3.2.** Suppose that Assumptions **B1** to **B5** or **B1** to **B4** and **B5'** hold true and $(Y_t)_{t \in \mathbb{Z}}$ is a strictly stationary process satisfying (5) and (6).

1) i) Assume that $\mathbb{E} (\|X_0\|_2^2) < +\infty$, then the estimators $\hat{\theta}_{\mathcal{P}\mathcal{L}}$, $\hat{\gamma}$ and $\hat{r}$ are strongly consistent, i.e.

$$ \hat{\theta}_{\mathcal{P}\mathcal{L}} \xrightarrow{T \to +\infty} \theta_0, \quad \hat{\gamma} \xrightarrow{T \to +\infty} \gamma_0 \quad \text{and} \quad \hat{r} \xrightarrow{T \to +\infty} r_0. $$

ii) Assume furthermore that $\gamma_0$ is located in the interior of $\Gamma$. Setting for $y_i \in \{0, 1\}$ and $s_i \in \mathbb{R}$, $h_{y_i}(s_i) = y_i \log \Phi(s_i) + (1 - y_i) \log \Phi(-s_i)$, we have the stochastic expansion

$$ \sqrt{T}(\hat{\gamma} - \gamma_0) = \frac{J_{\gamma_0}^{-1}}{\sqrt{T}} \sum_{t=p+1}^{T} \sum_{i=1}^{k} \dot{h}_{Y_{i,t}} (\lambda_{i,t}(\gamma_0)) \nabla \lambda_{i,t}(\gamma_0) + o_p(1), $$

where

$$ J_{\gamma_0} = - \sum_{i=1}^{k} \mathbb{E} [\ddot{h}_{Y_{i,0}} (\lambda_{i,0}(\gamma_0)) \nabla \lambda_{i,0}(\gamma_0) \nabla \lambda_{i,0}(\gamma_0)'] . $$

As a consequence, we have the asymptotic normality

$$ \sqrt{T}(\hat{\gamma} - \gamma_0) \xrightarrow{\mathcal{L}} \mathcal{N}_Q (0, J_{\gamma_0}^{-1} I_{\gamma_0} J_{\gamma_0}^{-1}), $$
where

$$I_{\gamma_0} = \sum_{i=1}^{k} \sum_{j=1}^{k} \mathbb{E} \left[ \hat{h}_{Y_{i,0}} (\lambda_{i,0}(\gamma_0)) \hat{h}_{Y_{j,0}} (\lambda_{j,0}(\gamma_0)) \nabla \lambda_{i,0}(\gamma_0) \nabla \lambda_{j,0}(\gamma_0)' \right].$$

2) Assume there exists $\kappa > 0$ such that for all $t \in \mathbb{Z}$

$$\mathbb{E} \left( \exp \left( \kappa \|X_t\|_2^2 \right) \right) < +\infty,$$

and assume furthermore that $\theta_0$ is located in the interior of $\Theta$.

i) Setting for any $t \in \mathbb{Z}$

$$m_t(\theta) = \log \left( \int_{\mathbb{R}^k} \prod_{i=1}^k \mathbb{1}_{Y_{i,t}} (\lambda_{i,t}(\gamma) + x_i) \varphi_R(x) dx \right),$$

we have the stochastic expansion

$$\sqrt{T} (\hat{\theta}_{\mathcal{PL}} - \theta_0) = \frac{J_{\theta_0}^{-1}}{\sqrt{T}} \nabla \mathcal{PL}(\theta_0) + o_{\mathbb{P}}(1)$$

where $J_{\theta_0} = -\mathbb{E}(\nabla^2 m_0(\theta_0))$.

As a consequence, we have the asymptotic normality

$$\sqrt{T} (\hat{\theta}_{\mathcal{PL}} - \theta_0) \xrightarrow{\mathcal{L}} \mathcal{N} (0, J_{\theta_0}^{-1} I_{\theta_0} J_{\theta_0}^{-1})$$

where $I_{\theta_0} = \mathbb{E}(\nabla m_0(\theta_0) \cdot \nabla m_0(\theta_0)')$.

ii) Setting

$$\mathbb{E}(\nabla^2 \ell_0(\theta_0)) = \mathbb{E} \left( \begin{array}{c|c} \frac{\partial^2 \ell_0}{\partial \gamma^2}(\theta_0) & \frac{\partial^2 \ell_0}{\partial \gamma \partial r}(\theta_0) \\ \hline \frac{\partial^2 \ell_0}{\partial r \partial \gamma}(\theta_0) & \frac{\partial^2 \ell_0}{\partial r^2}(\theta_0) \end{array} \right) = \left( \begin{array}{c|c} L_{1,1} & L_{2,1} \\ \hline L_{1,2} & L_{2,2} \end{array} \right)$$

we have the stochastic expansion

$$\sqrt{T}(\hat{r} - r_0) = -\frac{1}{\sqrt{T}} \sum_{t=p+1}^{T} L_{2,2}^{-1} \nabla_2 \ell_t(\theta_0) - \sqrt{T} L_{2,2}^{-1} L_{1,2} (\hat{\gamma} - \gamma_0) + o_{\mathbb{P}}(1).$$

In particular, $\sqrt{T}(\hat{\gamma} - \gamma_0, \hat{r} - r_0)$ also has a Gaussian limiting distribution with mean 0.

## Notes.

1. Our asymptotic results require some moment conditions on the exogenous covariates. For consistency, only a moment of order 2 is required. One can note that the covariates do not need to be bounded almost surely, which is important to be compatible with most of the time series models used in Econometrics. On the other hand, asymptotic normality of conditional likelihood estimators or of pairwise likelihood estimators of the correlation coefficients require an exponential moment for $\|X_0\|_2^2$. This condition is required to control the variance of the score which is a ratio of multiple integrals. Though restrictive, this condition is still satisfied for any stationary Gaussian time series. In both cases, bounded covariates sequences are allowed. An example of bounded time series in Ecology concerns the absence-presence of other species than the species of interest, as discussed in Section 5.2.
2. Although all the asymptotic variances may have a complicated expression and are not all defined explicitly in the previous theorem, all of them have a sandwich form. Note that confidence interval for model parameters can be also obtained by parametric bootstrap. The main difficulty for deriving asymptotic normality of the previous estimates is the presence of some parameters in the integration bounds of the multivariate Gaussian density, which leads to lengthy and tedious asymptotic expansions.

# 4 Parametric inference for panel type data

Our aim is now to generalize the results given in Section 3 to the case of panel data. Ecology is indeed a field where abundance data is collected from surveys at various locations and most of the time one or two times per year. Taking into account the information coming from different time points $t = 1, \dots, T$ and locations $j = 1, \dots, n$ is then primordial to get reasonable sample sizes.

Several time series $(X_{j,t}, Y_{j,t})_{1\le t\le T}$ are observed at different locations $j = 1, \dots, n$, with $X_{j,t}$ being a vector of exogenous covariates and $Y_{j,t}$ is the absence/presence vector observed at time $t$ and location $j$. We will assume that the same model (5), (6) is satisfied for all the observed time series. $M$-estimators are generally $\sqrt{nT}$-consistent in this context and the problem of statistical inaccuracy in our multivariate model due to short horizon times in many Ecological time series can be overcome if several dynamics are observed for a significant number of «individuals», which corresponds to ecosystems observed at different spatial locations in our context. Note that our models and results are by no means restricted to Ecological time series. In our real data example, the orders of magnitude of $n$ and $T$ are similar. Then for the asymptotics, we will work under a large $T$ and large $n$ scenario. Our joint limits will be obtained under the most general setup, that is when both $T$ and $n$ go to infinity but without other restrictions such as the standard diagonal path limit used in the panel data literature for which $T = T_n$ is such that $n/T \to c \neq 0$ as the index $n \to \infty$. We defer the reader to Hsiao (2022), Appendix 3B, for a description of various approaches for large $n$ and large $T$ asymptotics. We will however not follow the standard assumptions used for panels data with unobserved individual or time effects, though we will allow some observed random covariates to be only location-dependent or time-dependent. See Hahn and Kuersteiner (2011) or Fernández-Val and Weidner (2016) for asymptotic results for estimating regression parameters in panel data of this type in presence of unobserved heterogeneity. Additionally, in the latter references, weak law of large numbers for double-indexed partial sums are derived by combining mixing assumptions on the covariates and independence across individuals, conditionally on unobserved heterogeneity. Though general results for deriving mixing properties in discrete-valued time series models with covariates have been proposed recently in Truquet (2023), we will consider some strong (instead of weak) law of large numbers available without any stong mixing assumption.

## 4.1 Law of large numbers for partial sums of double-indexed sequences

We first present a rigorous version of the strong law of large numbers which is adapted to data collected at different time points $t$ and different locations $j$. Finding an appropriate result in the statistical literature is not straightforward since our aim here is to consider a complex modeling with covariates $X_{j,t}$ that can contain non-time varying components or some components not depending on the location $j$. Technically, the problem is to find an ergodic theorem adapted to multiparameter processes whose probability distributions are invariant under some specific transformations such as the time-shift.

There exists various formulations of the ergodic theorem for multiparameter processes. See for instance Nguyen and Zessin (1979) for a general ergodic theorem for spatial processes or Krengel (2011) for various ergodic theorems available in a general context. In this paper, we present here two versions of the strong law of large numbers which are based on an ergodic theorem given by Giap and Van Quang (2016), Lemma 3.1. A related result concerning more specifically stationary and ergodic random fields
can be also be found in Banys et al. (2010). An application of these results in the specific context of panel data is an important contribution of our work.

In what follows, we say that an array $(u_{n,m})_{n \in \mathbb{N}, m \in \mathbb{N}}$ converges toward $\ell$ when both $n$ and $m$ tend to $+\infty$ if

$$ \forall \varepsilon > 0, \exists k_0 \in \mathbb{N}, \forall n, m \geq k_0, |u_{n,m} - \ell| < \varepsilon. $$

We consider some strong laws of large numbers for an array $(Z_{j,t})_{j,t \geq 1}$ of the form

$$ Z_{j,t} = g(\zeta_{j,t}, \zeta_{j,t-1}, \dots), \quad \zeta_{j,t} = (X_{j,t-1}, \varepsilon_{j,t}) \qquad (11) $$

for some measurable mapping $g : (\mathbb{R}^d \times \mathbb{R}^k)^{\mathbb{N}} \to \mathbb{R}^{d'}$, $d'$ being a positive integer. Note that if $Z_{j,t} = (Y_{j,t}, Y_{j,t-1}, \dots, Y_{j,t-p})$ where $(Y_{j,t})_{t \in \mathbb{Z}}$, is solution of (5) and (6), such a representation is valid. From now on, for any $t \in \mathbb{Z}$ and $j \in \mathbb{N}^*$, we set

$$ \mathcal{F}_{j,t} = \sigma((X_{j,s}, \varepsilon_{j,s}) : s \leq t). \qquad (12) $$

Our first result is based on the following assumption.

**C1** For any $j \geq 1$, the stochastic process $(\zeta_{j,t})_{t \in \mathbb{Z}}$ is strictly stationary and its distribution does not depend on $j$. Moreover, setting $\mathcal{F}_{j,t} = \sigma((X_{j,s}, \varepsilon_{j,s}) : s \leq t)$ and $\mathcal{G}_j = \bigvee_{t \in \mathbb{Z}} \mathcal{F}_{j,t}$, the family of $\sigma$-fields $(\mathcal{G}_j)_{j \in \mathbb{N}^*}$ is mutually independent and for any $(j,t) \in \mathbb{N}^* \times \mathbb{Z}$, $\varepsilon_{j,t}$ is independent of $\mathcal{F}_{j,t-1}$.

**Proposition 4.1.** Suppose that (11) and C1 hold true. Then for all measurable mapping $f : \mathbb{R}^{d'} \to \mathbb{R}$ such that $\mathbb{E}(|f(Z_{1,1})| \times \log^+(|f(Z_{1,1})|)) < +\infty$, we have

$$ \lim_{n,T \to +\infty} \frac{1}{nT} \sum_{j=1}^{n} \sum_{t=1}^{T} f(Z_{j,t}) = \mathbb{E}(f(Z_{1,1})) \quad \text{a.s.} $$

### Notes.

1. Considering $\varepsilon_{j,t} = 0$ and setting $Z_{j,t} = X_{j,t}$, Proposition 4.1 shows that the strong law of large numbers for $(X_{j,t})_{j,t}$ is always valid when we have stationarity in the time direction and i.i.d. paths, up to an integrability condition which is slightly stronger than the usual one for stationary sequences. As discussed in Klesov and Molchanov (2017) in the independent setup, this moment condition is optimal and cannot be avoided.

2. Suppose that $X_{j,t} = (V_j, W_{j,t})$ takes values in $\mathbb{R}^{d_1} \times \mathbb{R}^{d_2}$ with $d = d_1 + d_2$, a useful scenario for considering both time-varying and non time-varying covariates. In the context of the absence/presence of waterbirds, Sebastian-Gonzalez et al. (2010) used ponds or vegetation characteristics for fixed covariates and absence/presence of other species across the time for time-varying covariates. Considering stationary but not necessarily ergodic paths $(Y_{j,t})_{t \in \mathbb{Z}}$ as in Theorem 2.1 is necessary if non-time varying covariates are included in the dynamic. To check C1, one possibility is to set $\bar{Z}_{j,t} = (W_{j,t}, \varepsilon_{j,t})$ and assume that the two families of random variables $\{V_j : j \geq 1\}$ and $\{\bar{Z}_{j,t} : j \geq 1, t \in \mathbb{Z}\}$ are independent, the $V_j$'s are i.i.d and $\{(\bar{Z}_{j,t})_{t \in \mathbb{Z}} : j \geq 1\}$ is a family composed of independent and identically distributed stationary stochastic processes. In a more subtle setup, the time-varying covariates $W_{j,t}$ may also depend on $V_j$. For instance, if $W_{j,t} = H(\eta_{j,t}, \eta_{j,t-1}, \dots)$ with $\eta_{j,t} = G(W_j, \kappa_{j,t})$ and $\{(\kappa_{j,t}, \varepsilon_{j,t}) : j \geq 1, t \in \mathbb{Z}\}$ is a family of i.i.d. random vectors independent from the sequence $(V_j)_{j \geq 1}$. In particular, the time-varying covariates $W_{j,t}$ can be represented as an infinite moving average sequence with innovation sequence $\eta_{j,t}$ that can still interact with $V_j$. One can then conclude that except the independence condition imposed for different $j$, the assumptions made on the covariates are sufficiently general to be compatible with practical applications.
3. Note that the integrability of $f \log^+(|f|)$ is guaranteed as soon as $|f|^{1+\delta}$ is integrable for some $\delta > 0$. For simplicity, this latter condition will be always assumed in the statements of our next results.

4. The strong law of large numbers might also be valid for some stationary random fields $(Z_{j,t})_{j \in \mathbb{N}, t \in \mathbb{Z}}$, without assuming independence between the sigma-fields $\mathcal{G}_j$. However, such results would require to specify the geometry of the domain indexing the spatial locations $j$. In panel data, independence assumption across the $j$ index is relevant when the different individuals do not interact. For spatial locations, this independence condition can be only considered as an approximation of the reality when the different locations are far enough from each other.

We also study another more specific scenario of the form (11) but without assuming independence across $j$. In this setup, the covariates have a general dynamical system representation with some common time-varying covariates.

**C'1** The covariates write as $X_{j,t} = (\Gamma_t, V_j, W_{j,t}) \in \mathbb{R}^{d_1} \times \mathbb{R}^{d_2} \times \mathbb{R}^{d_3}$, $d = d_1+d_2+d_3$. Moreover, there exist measurable mappings $H : (\mathbb{R}^{d_4})^N \to \mathbb{R}^{d_1}$, $L : (\mathbb{R}^{d_5})^N \to \mathbb{R}^{d_3}$ and $G : \mathbb{R}^{d_2} \times \mathbb{R}^{d_4} \times \mathbb{R}^{d_6} \to \mathbb{R}^{d_5}$ such that for $j \ge 1$ and $t \in \mathbb{Z}$,

$$ \Gamma_t = H(\eta_t, \eta_{t-1}, \dots), \quad W_{j,t} = L(\tilde{\eta}_{j,t}, \tilde{\eta}_{j,t-1}, \dots), \quad \tilde{\eta}_{j,t} = G(V_j, \eta_t, \kappa_{j,t}), \quad (13) $$

where $((\kappa_{j,t}, \varepsilon_{j,t}))_{j \in \mathbb{N}, t \in \mathbb{Z}}$, $(V_j)_{j \in \mathbb{N}}$ and $(\eta_t)_{t \in \mathbb{Z}}$ are three independent sequences of i.i.d. random variables taking values respectively in $\mathbb{R}^{d_6} \times \mathbb{R}^k$, $\mathbb{R}^{d_2}$ and $\mathbb{R}^{d_4}$.

**Proposition 4.2.** Suppose that (11) and C'1 hold true. Then for any measurable mapping $f : \mathbb{R}^{d'} \to \mathbb{R}$ such that

$$ \mathbb{E}(|f(Z_{1,1})| \cdot \log^+(|f(Z_{1,1})|)) < \infty, $$

we have

$$ \lim_{n,T \to +\infty} \frac{1}{nT} \sum_{j=1}^{n} \sum_{t=1}^{T} f(Z_{j,t}) = \mathbb{E}(f(Z_{1,1})) \quad \text{a.s.} $$

**Note.** It is not difficult to check that C'1 entails stationarity of the sequence $(\zeta_{j,t})_{t \in \mathbb{Z}}$, $j \ge 1$. This second scenario is interesting to take into account both local characteristics $(V_j, W_{j,t})$ of the system as well as global time-dependent and strictly exogenous characteristics $\Gamma_t$ (for instance, $\Gamma_t$ can represent the influence of macroeconomics variables in Econometrics or time-varying characteristics of the atmosphere in Ecology). Note that local covariates $W_{j,t}$ can be correlated with $\Gamma_t$ (through the impact of $\eta_t$) and can be sequentially exogenous and not necessarily strictly exogenous (since $\varepsilon_{j,t}$ and $\kappa_{j,t}$ are not necessarily independent). As assumed in (13), time-varying covariates can be represented as general random dynamical systems with correlated noise terms. Specification (13) includes general moving averages and many autoregressive processes. For instance $\Gamma_t = M\Gamma_{t-1} + \eta_t$, $W_{j,t} = NW_{j,t-1} + \tilde{\eta}_{j,t}$, where $M$ and $N$ are square matrices of appropriate dimension and spectral radius less than 1.

## 4.2 Asymptotic results for panel data

Before generalizing our estimation results from section 3, let us specify the framework of our statistical model for panel data.

In this entire section, we consider for $j \in \mathbb{N}^*$ two sequences $(X_{j,t})_{t \in \mathbb{Z}}$ and $(\varepsilon_{j,t})_{t \in \mathbb{Z}}$ of random vectors taking values in $\mathbb{R}^d$ and $\mathbb{R}^k$ respectively. We then define an array of random binary vectors $(Y_{j,t})_{j \in \mathbb{N}^*, t \in \mathbb{Z}}$ such that for any $j \in \mathbb{N}^*, t \in \mathbb{Z}$, and $i \in \{1, \dots, k\}$

$$ Y_{i,j,t} = \mathbb{1}_{(0, +\infty)} (\lambda_{i,j,t} + \varepsilon_{i,j,t}), \quad (14) $$
with

$$ \lambda_{j,t} = C + \sum_{l=1}^{p} A_l Y_{j,t-l} + B X_{j,t-1}. \quad (15) $$

**Note.** All matrices $A_1, \dots, A_p$ belong to $\mathbb{R}^{k \times k}$, $B \in \mathbb{R}^{k \times d}$, and $C \in \mathbb{R}^k$. These matrices are fixed in equation (15) and do not depend on $j$.

In the literature of panel data, it is classical to have heterogeneity across individual $j$ by adding an intercept $c_j$ to equation (15). The estimation of parameters $A_1, \dots, A_p$ and $B$ is then the central problem. See for instance Hsiao (2022), Chapter 7, for a classical reference on this topic. However, in the context of discrete data, there is no generic way for consistently estimating the fixed regression parameters. For instance, considering $(Y_{j,t})_{1 \le j \le N, 1 \le t \le T}$ for small $T$ and large $N$ in univariate binary models, the MLE of regression parameters can be inconsistent and an alternative consistent estimator is known in the logistic case, but not in the probit case. In this paper, we do not investigate such unobserved heterogeneity problem, but we allow observable random covariates which can vary only across individuals or only across time and which can play the rule of an observed heterogeneity.

We consider the following assumptions which are the analogs of the single path analysis.

**C2** For any $j \ge 1$ and $t \in \mathbb{Z}$, the random vector $\varepsilon_{j,t}$ are normally distributed

$$ \varepsilon_{j,t} \sim \mathcal{N}(0, R) $$

where $R$ is a correlation matrix, i.e.

$$ \forall i_1, i_2 \in \{1, \dots, k\}, \quad -1 \le R(i_1, i_2) \le 1 $$

and $R(i, i) = 1$ for any $i \in \{1, \dots, k\}$.

**C3** All matrices $A_1, \dots, A_p$ belong to a compact set $\mathcal{A} \subset \mathbb{R}^{k \times k}$, $B$ belongs to a compact set $\mathcal{B} \subset \mathbb{R}^{k \times d}$, $C$ belongs to a compact subset of $\mathbb{R}^k$ and $R$ belongs to a compact set $\mathcal{R}$ included in the set of positive definite symmetric matrices in $\mathbb{R}^{k \times k}$.

**C4** For any couple $(j, t) \in \mathbb{N}^* \times \mathbb{Z}$ and all $v \in \mathbb{R}^d$, we have for any $\sigma(\varepsilon_{j,t}) \lor \mathcal{F}_{j,t-1}$-measurable application $G$

$$ \mathbb{P}(v' \cdot X_{j,t} = G) < 1. $$

Keeping the same notations of Section 3.2, we finally focus on the estimation of the true vector of parameters $\theta_0 = (\gamma_0, r'_0)$ satisfying (14) and (15). We rewrite the two estimation procedures mentioned in section 3.2 in the specific framework of panel data.

Note that by assumption **C3**, $\theta$ belongs to a compact set $\Theta$, $\gamma$ belongs to a compact set $\Gamma$, and without ambiguity in the notations, $r$ belongs to a compact set $\mathcal{R} = \prod_{1 \le i_1 < i_2 \le k} \mathcal{R}_{i_1, i_2}$, where each set $\mathcal{R}_{i_1, i_2}$ is a compact set included in $(-1, 1)$.

1) We can first consider the so-called pseudo conditional log-likelihood. For any $t \in \mathbb{Z}$ and $n \in \mathbb{N}^*$, the distribution of $(Y_{1,t}, \dots, Y_{n,t})$ conditionally to $\mathcal{F}_{t-1}$ is given by

$$ \mathcal{L}((Y_{1,t}, \dots, Y_{n,t}) | \bigvee_{j=0}^{n} \mathcal{F}_{j,t-1}) = \sum_{s \in \{0,1\}^n} \left\{ \prod_{j=1}^{n} \left( \int_{\mathbb{R}^k} \prod_{i=1}^{k} 1_{I_{Y_{i,j,t}}}(\lambda_{i,j,t}(\gamma) + x_i) \varphi_R(x) dx \right) \right\} \delta_s. $$

Similarly to the single path case, we can define the pseudo conditional log-likelihood function by

$$ \mathcal{PL}(\theta) = \sum_{t=p+1}^{T} \sum_{j=1}^{n} \log \left( \int_{\mathbb{R}^k} \prod_{i=1}^{k} 1_{I_{Y_{i,j,t}}}(\lambda_{i,j,t}(\gamma) + x_i) \varphi_R(x) dx \right), \quad (16) $$
and we set $\hat{\theta}_{\mathcal{PL}} = \underset{\theta \in \Theta}{\operatorname{argmax}} \mathcal{PL}(\theta)$.

Once again, one can note that (16) is based on the density of $((Y_{1,t}, \dots, Y_{n,t}))_{p+1 \le t \le T}$ conditionally on $(Y_{1,1}, \dots, Y_{n,1}), \dots, (Y_{1,p}, \dots, Y_{n,p})$ and $(X_{1,1}, \dots, X_{n,1}), \dots, (X_{1,T}, \dots, X_{n,T})$, but such a density coincides with (16) only if the $X_{j,t}$'s are independent from the $\varepsilon_{j,t}$'s, which is not necessarily the case in our framework.

2) It is also possible to consider a two-step method as in Section 3.2, by first estimating the vector of parameters $\gamma$, then the correlation matrix $R$. We thus define

$$ \begin{aligned} \overline{\mathcal{PL}}(\gamma) &= \sum_{t=p+1}^{T} \sum_{j=1}^{n} \sum_{i=1}^{k} (Y_{i,j,t} \log [1 - \Phi(-\lambda_{i,j,t}(\gamma))] + (1 - Y_{i,j,t}) \log [\Phi(-\lambda_{i,j,t}(\gamma))]) \\ &= \sum_{t=p+1}^{T} \sum_{j=1}^{n} \sum_{i=1}^{k} (Y_{i,j,t} \log [\Phi(\lambda_{i,j,t}(\gamma))] + (1 - Y_{i,j,t}) \log [\Phi(-\lambda_{i,j,t}(\gamma))]) \end{aligned} $$

and we set $\hat{\gamma} = \underset{\gamma \in \Gamma}{\operatorname{argmax}} \overline{\mathcal{PL}}(\gamma)$.

In order to obtain an estimate for $R$, it is still possible to plug the estimate $\hat{\gamma}$ of $\gamma$ in the expression of function (16), and maximize it with respect to $R$, but it is still complicated numerically speaking.

A pairwise conditional likelihood approach similar to the one in section 3.2 is also possible, and setting for $i_1 < i_2$, $R_{i_1, i_2} = \begin{pmatrix} 1 & r_{i_1, 2} \\ r_{i_1, i_2} & 1 \end{pmatrix}$ where $r_{i_1, i_2} = R(i_1, i_2)$, we can define

$$ \begin{aligned} \hat{r}_{i_1, i_2} &= \underset{r_{i_1, i_2} \in \mathcal{R}_{i_1, i_2}}{\operatorname{argmax}} \sum_{t=p+1}^{T} \sum_{j=1}^{n} \log \left( \int_{I_{Y_{i_1, j,t}} - \lambda_{i_1, j,t}(\hat{\gamma})} \int_{I_{Y_{i_2, j,t}} - \lambda_{i_2, j,t}(\hat{\gamma})} \varphi_{R_{i_1, i_2}}(x_{i_1}, x_{i_2}) dx_{i_1} dx_{i_2} \right) \\ &= \underset{r_{i_1, i_2} \in \mathcal{R}_{i_1, i_2}}{\operatorname{argmax}} \sum_{t=p+1}^{T} \sum_{j=1}^{n} \log \left( \int_{I_{Y_{i_1, j,t}} - \lambda_{i_1, j,t}(\hat{\gamma})} \Phi \left( (2Y_{i_2, j,t} - 1) \frac{(\lambda_{i_2, j,t}(\hat{\gamma}) + r_{i_1, i_2} x_{i_1})}{\sqrt{1 - r_{i_1, i_2}^2}} \right) \varphi(x_{i_1}) dx_{i_1} \right). \end{aligned} $$

Note that once again, under assumption C4, for any $t \in \mathbb{Z}$, any $n \in \mathbb{N}^*$ and $s_{i_1, 1}, \dots, s_{i_1, n}, s_{i_2, 1}, \dots, s_{i_2, n} \in \{0, 1\}$,

$$ \begin{aligned} & \prod_{j=1}^{n} \int_{I_{s_{i_1,j}-\lambda_{i_1,j,t}(\gamma_0)}} \Phi \left( (2s_{i_2,j} - 1) \frac{\lambda_{i_2,j,t}(\gamma) + r_{i_1,i_2}e_{i_1}}{\sqrt{1-r_{i_1,i_2}^2}} \right) \varphi(e_{i_1})de_{i_1} \\ &= \mathbb{P}((Y_{i_1,1,t}, Y_{i_2,1,t}) = (s_{i_1,1}, s_{i_2,1}), \dots, (Y_{i_1,n,t}, Y_{i_2,n,t}) = (s_{i_1,n}, s_{i_2,n}) | \mathcal{F}_{t-1}), \end{aligned} $$

which explains the terminology pairwise conditional likelihood, even it still doesn't correspond to a proper likelihood function.

Similarly as in the single path case, we denote $\hat{r} = (\hat{r}_{1,2}, \dots, \hat{r}_{1,k}, \hat{r}_{2,3}, \dots, \hat{r}_{2,k}, \dots \hat{r}_{k-1,k})$, and we set for any $j \in \mathbb{N}^*, t \in \mathbb{Z}$ and $\theta = (\gamma, r') \in \Theta$

$$ \ell_{j,t}(\theta) = \sum_{1 \le i_1 < i_2 \le k} \log \left( \int_{I_{Y_{i_1,j,t}} - \lambda_{i_1,j,t}(\gamma)} \Phi \left( (2Y_{i_2,j,t} - 1) \frac{\lambda_{i_2,j,t}(\gamma) + r_{i_1,i_2}x_{i_1}}{\sqrt{1 - r_{i_1,i_2}^2}} \right) \varphi(x_{i_1}) dx_{i_1} \right). $$

**Theorem 4.3.** Suppose that either Assumption C1 or Assumption C'1 is valid and Assumptions C2-C3-C4 hold true. For any $j \in \{1, \dots, n\}$, there then exists a unique strictly stationary process $(Y_{j,t})_{t \in \mathbb{Z}}$ satisfying (14) and (15). Moreover, the following assertions are valid.
1) i) Assume that for all $j \in \mathbb{N}^*$ and $t \in \mathbb{Z}$, $\mathbb{E} (\|X_{j,t}\|_2^{2+\delta}) < +\infty$ for some $\delta > 0$, then the estimators $\hat{\theta}_{\mathcal{P}\mathcal{L}}$, $\hat{\gamma}$ and $\hat{r}$ are strongly consistent, i.e.

$$ \hat{\theta}_{\mathcal{P}\mathcal{L}} \underset{n,T \to \infty}{\longrightarrow} \theta_0, \quad \hat{\gamma} \underset{n,T \to \infty}{\longrightarrow} \gamma_0 \quad \text{and} \quad \hat{r} \underset{n,T \to \infty}{\longrightarrow} r_0. $$

ii) Assume furthermore that $\gamma_0$ is located in the interior of $\Gamma$. Setting for $y_i \in \{0,1\}$ and $s_i \in \mathbb{R}$, $h_{y_i}(s_i) = y_i \log(\Phi(s_i)) + (1 - y_i) \log(\Phi(-s_i))$, we have the stochastic expansion

$$ \sqrt{nT}(\hat{\gamma} - \gamma_0) = \frac{J_{\gamma_0}^{-1}}{\sqrt{nT}} \sum_{j=1}^{n} \sum_{t=p+1}^{T} \hat{h}_{Y_{i,j,t}}(\lambda_{i,j,t}(\gamma_0)) \nabla \lambda_{i,j,t}(\gamma_0) + o_{\mathbb{P}}(1), $$

when $n, T \to \infty$, and where

$$ J_{\gamma_0} = - \sum_{i=1}^{k} \mathbb{E} \left[ \ddot{h}_{Y_{i,1,1}}(\lambda_{i,1,1}(\gamma_0)) \nabla \lambda_{i,1,1}(\gamma_0) \nabla \lambda_{i,1,1}(\gamma_0)' \right]. $$

As a consequence, we have the asymptotic normality

$$ \sqrt{nT}(\hat{\gamma} - \gamma_0) \xrightarrow{\mathcal{L}} \mathcal{N}(0, J_{\gamma_0}^{-1} I_{\gamma_0} J_{\gamma_0}^{-1}) $$

when $n, T \to \infty$, where

$$ I_{\gamma_0} = \sum_{i=1}^{k} \sum_{j=1}^{k} \mathbb{E} \left[ \hat{h}_{Y_{i,1,1}}(\lambda_{i,1,1}(\gamma_0)) \hat{h}_{Y_{j,1,1}}(\lambda_{j,1,1}(\gamma_0)) \nabla \lambda_{i,1,1}(\gamma_0) \nabla \lambda_{j,1,1}(\gamma_0)' \right]. $$

2) Assume there exists $\kappa > 0$ such that for all $j \in \mathbb{N}^*$ and $t \in \mathbb{Z}$

$$ \mathbb{E} (\exp (\kappa \|X_{j,t}\|_2^2)) < +\infty, \qquad (17) $$

and assume furthermore that $\theta_0$ is located in the interior of $\Theta$.

i) Setting for any $j \in \mathbb{N}^*$ and $t \in \mathbb{Z}$

$$ m_{j,t}(\theta) = \log \left( \int_{\mathbb{R}^k} \prod_{i=1}^k \mathbb{1}_{Y_{i,j,t}} (\lambda_{i,j,t}(\gamma) + x_i) \varphi_R(x) dx \right), $$

we have the stochastic expansion

$$ \sqrt{nT} (\hat{\theta}_{\mathcal{P}\mathcal{L}} - \theta_0) = \frac{J_{\theta_0}^{-1}}{\sqrt{nT}} \nabla \mathcal{P}\mathcal{L}(\theta_0) + o_{\mathbb{P}}(1) $$

when $n, T \to +\infty$, and where $J_{\theta_0} = -\mathbb{E}(\nabla^2 m_{1,1}(\theta_0))$.

As a consequence, we have the asymptotic normality

$$ \sqrt{nT} (\hat{\theta}_{\mathcal{P}\mathcal{L}} - \theta_0) \xrightarrow{\mathcal{L}} \mathcal{N} (0, J_{\theta_0}^{-1} I_{\theta_0} J_{\theta_0}^{-1}) $$

when $n, T \to +\infty$, where $I_{\theta_0} = \mathbb{E}(\nabla m_{1,1}(\theta_0) \cdot \nabla m_{1,1}(\theta_0))'$.

ii) Setting

$$ \mathbb{E}(\nabla^2 \ell_0(\theta_0)) = \mathbb{E} \left( \begin{array}{c|c} \frac{\partial^2 \ell_0}{\partial \gamma^2}(\theta_0) & \frac{\partial^2 \ell_0}{\partial \gamma \partial r}(\theta_0) \\ \hline \frac{\partial^2 \ell_0}{\partial r \partial \gamma}(\theta_0) & \frac{\partial^2 \ell_0}{\partial r^2}(\theta_0) \end{array} \right) = \left( \begin{array}{c|c} L_{1,1} & L_{2,1} \\ \hline L_{1,2} & L_{2,2} \end{array} \right) \qquad (18) $$
we have the stochastic expansion

$$ \sqrt{nT}(\hat{r} - r_0) = -\frac{1}{\sqrt{nT}} \sum_{t=p+1}^{T} \sum_{j=1}^{n} L_{2,2}^{-1} \nabla_2 \ell_{j,t}(\theta_0) - \sqrt{nT} L_{2,2}^{-1} L_{1,2}(\hat{\gamma} - \gamma_0) + o_{\mathbb{P}}(1) $$

when $n, T \to +\infty$.

In particular, $\sqrt{nT}(\hat{\gamma} - \gamma_0, \hat{r} - r_0)$ also has a Gaussian limiting distribution with mean 0.

### Notes.

1. As expected, the convergence rate of the different estimators is $\sqrt{nT}$ and $(n, T)$ can grow to infinity with some arbitrary rates.

2. The law of large numbers derived for two-parameters partial sums plays a central role in the proof of Theorem 4.3. With respect to the single path analysis, we simply need to reinforce the assumption of square-integrability of the covariates by the existence of a moment greater than 2. In the supplementary file, we provide some general asymptotic results for parametric M-estimators, in the setup of longitudinal observations and under Assumption C1 or Assumption C'1. These generic results can be also applied to general autoregressive time series models with exogenous covariates.

# 5 Numerical applications

In this section, we first compare the two inference strategies mentioned in the previous section, in the case of simulated data. We will then apply the two-steps method to a real dataset of fish catches in the North Sea, in order to get some insights about the interactions between some particular species.

## 5.1 Simulated data and methods comparison

Our aim here is to compare the performances of our two estimation strategies. To do so, we perform 1000 simulations of an absence-presence process of two species, depending on one exogenous variable. For each simulation, the time dimension is fixed at $T = 100$, and the number of time series paths is $n = 50$. For simulation number $s \in \{1, \dots, 1000\}$, we thus denote $(Y_{j,t}^{(s)})_{t \in \mathbb{Z}}$ the time series valued in $\{0, 1\}^2$ for path $j \in \{1, \dots, 50\}$.

Furthermore, we denote $(X_{j,t}^{(s)})_{t \in \mathbb{Z}}$ the process of covariates valued in $\mathbb{R}$ that impacts the process $(Y_{j,t}^{(s)})_{t \in \mathbb{Z}}$ in our study. For any $j \in \{1, \dots, 50\}$, $X_{j,t}^{(s)}$ is a stationary ARMA(3, 1) process.

The covariates process is centered and standardized before simulating the process $(Y_{j,t}^{(s)})_{1 \le t \le T}$ according to the following equation

$$ \forall t \in \mathbb{Z}, \forall i \in \{1, 2\}, \quad Y_{i,j,t}^{(s)} = 1_{(0, +\infty)} \left( \lambda_{i,j,t}^{(s)} + \varepsilon_{i,j,t}^{(s)} \right) $$

where $\lambda_{j,t}^{(s)} = A \cdot Y_{j,t-1}^{(s)} + B \cdot X_{j,t-1}^{(s)} + C$ with

$$ A = \begin{pmatrix} 0.3 & -0.5 \\ 0.2 & 0.7 \end{pmatrix}, \quad B = \begin{pmatrix} -0.5 \\ 0.6 \end{pmatrix}, \quad C = \begin{pmatrix} 0.2 \\ 0.4 \end{pmatrix} $$

and the variables $\varepsilon_{j,t}^{(s)}$ are independent and identically distributed, $\varepsilon_{j,t}^{(s)} \sim \mathcal{N}(0, R)$, where the correlation matrix $R$ is given by

$$ R = \begin{pmatrix} 1 & -0.2 \\ -0.2 & 1 \end{pmatrix}. $$
One can find in Tables 1 and 2 our estimation results with both 1-step and 2-steps methods discussed in Section 4.

For the initialization of coefficient $A(k, l)$ in the optimization procedure, we actually compute over all 50 paths how many times the presence of species $l$ at time $t$ is followed by the presence of species $k$ at time $t + 1$, and how many times it is not. These two latter sums are pondered by 1 or -1 whether we consider the positive or negative impact of species $l$ on species $l$. We then divide the obtained result by the number of time steps over all 50 sites to get an initial guess for coefficient $A(k, l)$.

The initialization of coefficient $B(k)$ is simply given by the mean over the 50 paths observed of the correlation estimated between $Y_{k,j,t}^{(s)}$ and $X_{j,t-1}^{(s)}$. For each estimation of the correlation, we proceed here as if for $1 \le k \le 50$, $(Y_{k,j,t}^{(s)}, X_{j,t}^{(s)})_{1 \le t \le 100}$ were i.i.d. random variables.

Coefficients of type $C(k)$ are simply initialized by the proportion $p_k$ of presence over time $t$ and space $j$ for species $k$, which is then transformed by $2 \times p_k - 1$.

Finally, for a coefficient $R(k, l)$, we simply estimate the correlation between $Y_{k,j,t}^{(s)}$ and $Y_{l,j,t}^{(s)}$. Here, we also proceed as if the random variables $(Y_{k,j,t}^{(s)}, Y_{l,j,t}^{(s)})_{1 \le t \le 100 \atop 1 \le j \le 50}$ were i.i.d.

In order to compare the two methods, we fixed in both cases the same stopping criterion to end the numeric optimization, and used the same 100 CPUs. The estimates of $\hat{\gamma}$ in the 2-steps method were obtained using the Nelder-Mead simplex algorithm (see Nelder and Mead (1965)), while the estimates of $\hat{R}$ in the 2-steps method, as well as the estimates of $\hat{\theta}$ in the 1-step method, were computed with the COBYLA algorithm (see Powell (1994)) due to the constraints on the correlation matrix.

Table 1 and Table 2 exhibit the bias and variance of the different estimates. Note that for this moderate sample size, the 2-steps method performed significantly better both in terms of bias and MSE, while the variances of the estimates are very similar for the two methods. We think that this larger bias can be due to numerical issues. Indeed, the 1-step method is very difficult to implement with a numerical approximation of multiple integrals. Moreover, the presence of many local maxima makes the optimization of (16) quite difficult. The computation time for the 2-steps method was about 21 minutes, while it was about 2 hours and 30 minutes for the 1-step method. Such simulation shows that the 2-steps method is an interesting choice since it reduces the computational cost of the 1-step method and is still accurate for moderate sample sizes. Only the 2-step method will be used for the real data application given in the next section.

<table><caption>Table 1: Estimates Results for the 1-step Method. The mean is taken over the 1000 simulations</caption><thead><tr><th>Parameter</th><th>True Value</th><th>Mean Estimate</th><th>MSE</th><th>Bias</th><th>Variance</th></tr></thead><tbody><tr><td>A(1, 1)</td><td>0.3</td><td>0.358</td><td>0.005</td><td>0.058</td><td>0.002</td></tr><tr><td>A(1, 2)</td><td>-0.5</td><td>-0.338</td><td>0.283</td><td>0.162</td><td>0.002</td></tr><tr><td>A(2, 1)</td><td>0.2</td><td>0.186</td><td>0.002</td><td>-0.014</td><td>0.002</td></tr><tr><td>A(2, 2)</td><td>0.7</td><td>0.680</td><td>0.003</td><td>-0.020</td><td>0.003</td></tr><tr><td>B(1)</td><td>-0.5</td><td>-0.411</td><td>0.001</td><td>-0.011</td><td>0.001</td></tr><tr><td>B(2)</td><td>0.6</td><td>0.601</td><td>0.001</td><td>0.001</td><td>0.001</td></tr><tr><td>C(1)</td><td>0.2</td><td>0.014</td><td>0.036</td><td>-0.186</td><td>0.001</td></tr><tr><td>C(2)</td><td>0.4</td><td>0.427</td><td>0.004</td><td>0.027</td><td>0.003</td></tr><tr><td>R(1, 2)</td><td>-0.2</td><td>-0.199</td><td>0.001</td><td>-0.001</td><td>0.001</td></tr></tbody></table>
<table><caption>Table 2: Estimates Results for the 2-steps Method. The mean is taken over the 1000 simulations</caption><thead><tr><th>Parameter</th><th>True Value</th><th>Mean Estimate</th><th>MSE</th><th>Bias</th><th>Variance</th></tr></thead><tbody><tr><td>A(1, 1)</td><td>0.3</td><td>0.298</td><td>0.002</td><td>-0.002</td><td>0.002</td></tr><tr><td>A(1, 2)</td><td>-0.5</td><td>-0.503</td><td>0.003</td><td>-0.003</td><td>0.003</td></tr><tr><td>A(2, 1)</td><td>0.2</td><td>0.200</td><td>0.002</td><td>0.000</td><td>0.002</td></tr><tr><td>A(2, 2)</td><td>0.7</td><td>0.700</td><td>0.003</td><td>0.000</td><td>0.003</td></tr><tr><td>B(1)</td><td>-0.5</td><td>-0.400</td><td>0.001</td><td>0.100</td><td>0.001</td></tr><tr><td>B(2)</td><td>0.6</td><td>0.600</td><td>0.001</td><td>0.000</td><td>0.001</td></tr><tr><td>C(1)</td><td>0.2</td><td>0.203</td><td>0.003</td><td>0.003</td><td>0.003</td></tr><tr><td>C(2)</td><td>0.4</td><td>0.401</td><td>0.003</td><td>0.001</td><td>0.003</td></tr><tr><td>R(1, 2)</td><td>-0.2</td><td>-0.199</td><td>0.001</td><td>0.001</td><td>0.001</td></tr></tbody></table>

## 5.2 Fishing data

The following example relies on a dataset of fish catches provided publicly by the website of the International Council for the Exploration of the Sea (ICES), Copenhagen. More precisely, this dataset contains the Catch Per Unit Effort (CPUE) per length per subarea from the North Sea International Bottom Trawl Survey (NS-IBTS), and was downloaded the 29th may 2024 (ICES, 2024).

The CPUE corresponds here to the number of catches per hour, and is somehow an indicator of abundance for each fish species.

The downloaded NS-IBTS dataset thus contains the CPUE observed for each subarea in the North Sea (see Figure 1) in Area 1, from 1971 to 2024. There are nine fishery species studied in this dataset, scientifically named *Clupea Harengus*, *Gadus Morhua*, *Melanogrammus Aeglefinus*, *Merlangius Merlangus*, *Pleuronectes Platessa*, *Pollachius Virens*, *Scomber Scombrus*, *Sprattus Sprattus* and *Trisopterus Esmarkii*, which are subdivided in different length classes.

<figure><img src="image_2.png" /><figcaption>Figure 1: Areas and subareas of the NS-IBTS dataset. The dots correspond to the 36 subareas used in our study.</figcaption></figure>
For practical reasons, our study focuses on the time interval 1975-2024, and we first aggregate the CPUEs allocated by length for each species, then we remove the subareas containing too much missing data, which leaves us with 36 subareas (see Figure 1), corresponding to 36 paths over time for our absence-presence process. Moreover, we compute the first tiercile of the CPUEs for each species over time and space, in order to set the absence-presence index to 0 if the CPUE observed is smaller than this tiercile, and we set the absence-presence index to 1 otherwise. Hence, a value 0 does not necessarily mean that the species is absent, but that the species is somehow less abundant than usual. Finally, we impute the few remaining missing data at random, according to a Bernoulli distribution $B(p)$ where $p$ is the proportion of ones obtained in the whole dataset for the species for which we want to impute the absence-presence.

We decide to focus on the following three particular species: Gadus Morhua, Merlangius Merlangus and Sprattus Sprattus. The processes of absence-presence for the six other species are used here as exogenous variables. The model we fit to our data is thus

$$ \forall t \in \mathbb{Z}, \forall j \in \{1, \dots, 36\}, \forall i \in \{1, 2, 3\}, Y_{i,j,t} = 1_{(0, +\infty)} (\lambda_{i,j,t} + \varepsilon_{i,j,t}), $$

where $\lambda_{j,t} = A \cdot Y_{j,t-1} + B \cdot X_{j,t-1} + C$ with $A \in \mathbb{R}^{3 \times 3}$, $B \in \mathbb{R}^{3 \times 6}$, $X_{j,t-1} \in \{0, 1\}^6$ and $C \in \mathbb{R}^3$. $X_{j,t}$ corresponds to the vector of covariates, i.e. the absence-presence process of the six non-studied species, and the $\varepsilon_{j,t}$ are i.i.d with normal distribution $N(0, R)$, $R \in \mathbb{R}^{3 \times 3}$ being a correlation matrix. By doing so, we actually obtain a sub-model of the complete model with nine species, but with no covariates.

Due to the long time required by the 1-step estimation strategy, we only estimate the values of parameters $A, B, C$ and $R$ with the 2-steps method. Estimation results are stored in Table 3.

Since we are interested in the interpretation of the model's parameters, we then want to test the significance of the estimated parameters. Denoting

$$ \theta = (A(1, 1), \dots, A(3, 3), B(1, 1), \dots, B(3, 6), C(1), C(2), C(3), R(2, 1), R(3, 1), R(3, 2))', $$

$\theta_0$ the true value of the parameter $\theta$ and $\hat{\theta}$ its estimation, we have by Theorem 4.3

$$ \sqrt{nT} (\hat{\theta} - \theta_0) \xrightarrow{\mathcal{L}} \mathcal{N}(0, V_{\theta_0}) $$

where $V_{\theta_0}$ is an unknown covariance matrix. However, the quantiles of the distribution $N(0, V_{\theta_0})$ can be approximated by a bootstrap estimation (see DiCiccio and Efron (1996) for a review of bootstrap confidence intervals). To do so, we simulate a thousand absence-presence processes with 36 paths and duration 50 with the estimated parameter $\hat{\theta}$. These simulations then give us a thousand estimates of $\hat{\theta}$, denoted $\hat{\theta}_1^*, \dots, \hat{\theta}_{1000}^*$.

Computing the empirical quantiles of these latter estimates, we obtain good approximations of the quantiles of the limit distribution $N(0, V_{\theta_0})$.

We can then deduce easily some confidence intervals for the coefficients of parameter $\theta_0$. These results are given in Table 3.

Finally, we propose some interpretations concerning the results obtained. For example, coefficient $A(3, 2) \simeq 0.393$ being significantly positive, it suggests that the presence of species 2 (Merlangius merlangus) at time $t-1$ has a positive impact on the presence of species 2 (Sprattus sprattus) at time $t$. Indeed, it appears in a survey (see CEFAS (2024), consulted online the 10th of June, 2024) conducted by the Center for Environment Fishes and Aquaculture Science (CEFAS), monitoring fishes stomachs, that species Sprattus Sprattus is a particular predator for the eggs of species like Merlangius merlangus.

Coefficient $A(2, 3) \simeq 0.448$ is however significantly positive, meaning that the presence of species Sprattus sprattus at time $t-1$ has also a positive impact on the presence of species Merlangius merlangus at time $t$, which might seem in contradiction with the previous interpretation. This actually
discloses an underlying phenomenon that allows these two species to live in symbiosis. This insight is indeed confirmed by the estimated coefficient $R(2, 3) \simeq 0.158$, which is also significantly positive.

The model proposed in this paper also allows to give some insights about fishes interactions that might not be obvious at first sight. This is the case for species *Scomber scombrus* and *Merlangius merlangus*, as coefficient $B(2, 5) \simeq -0.184$ is significantly negative. This means that *Scomber scombrus* has a negative impact on the presence of *Merlangius merlangus*, although there is no direct indication that these two predator species are in competition.

Many other interpretations could be deduced from Table 3, but this is actually out of the scope of this paper, where the numerical study is just here given as an illustration.

<table><caption>Table 3: Estimation results for the NS-IBTS dataset.</caption><thead><tr><th>Parameter</th><th>Estimate</th><th>95% Confidence Interval</th></tr></thead><tbody><tr><td>A(1, 1)</td><td>0.841</td><td>[0.715;0.982]</td></tr><tr><td>A(1, 2)</td><td>0.313</td><td>[0.173;0.443]</td></tr><tr><td>A(1, 3)</td><td>0.218</td><td>[-0.017;0.430]</td></tr><tr><td>A(2, 1)</td><td>0.144</td><td>[0.005;0.290]</td></tr><tr><td>A(2, 2)</td><td>1.004</td><td>[0.861;1.135]</td></tr><tr><td>A(2, 3)</td><td>0.448</td><td>[0.197;0.679]</td></tr><tr><td>A(3, 1)</td><td>0.117</td><td>[-0.062;0.297]</td></tr><tr><td>A(3, 2)</td><td>0.393</td><td>[0.195;0.576]</td></tr><tr><td>A(3, 3)</td><td>0.000</td><td>[-0.002;0.002]</td></tr><tr><td>B(1, 1)</td><td>-0.066</td><td>[-0.201;0.079]</td></tr><tr><td>B(1, 2)</td><td>0.122</td><td>[-0.016;0.272]</td></tr><tr><td>B(1, 3)</td><td>0.065</td><td>[-0.074;0.205]</td></tr><tr><td>B(1, 4)</td><td>0.250</td><td>[0.086;0.401]</td></tr><tr><td>B(1, 5)</td><td>0.060</td><td>[-0.094;0.223]</td></tr><tr><td>B(1, 6)</td><td>0.103</td><td>[-0.043;0.266]</td></tr><tr><td>B(2, 1)</td><td>0.513</td><td>[0.343;0.656]</td></tr><tr><td>B(2, 2)</td><td>0.233</td><td>[0.091;0.379]</td></tr><tr><td>B(2, 3)</td><td>0.198</td><td>[0.046;0.334]</td></tr><tr><td>B(2, 4)</td><td>-0.198</td><td>[-0.352;-0.041]</td></tr><tr><td>B(2, 5)</td><td>-0.184</td><td>[-0.355;-0.003]</td></tr><tr><td>B(2, 6)</td><td>0.094</td><td>[-0.052;0.241]</td></tr><tr><td>B(3, 1)</td><td>0.004</td><td>[-0.025;0.027]</td></tr><tr><td>B(3, 2)</td><td>0.023</td><td>[-0.122;0.152]</td></tr><tr><td>B(3, 3)</td><td>-0.022</td><td>[-0.166;0.119]</td></tr><tr><td>B(3, 4)</td><td>-0.009</td><td>[-0.053;0.056]</td></tr><tr><td>B(3, 5)</td><td>-0.017</td><td>[-0.102;0.091]</td></tr><tr><td>B(3, 6)</td><td>0.093</td><td>[-0.122;0.250]</td></tr><tr><td>C(1)</td><td>-0.644</td><td>[-0.823;-0.452]</td></tr><tr><td>C(2)</td><td>-0.708</td><td>[-0.873;-0.505]</td></tr><tr><td>C(3)</td><td>-1.697</td><td>[-1.889;-1.442]</td></tr><tr><td>R(1, 2)</td><td>0.229</td><td>[0.200;0.407]</td></tr><tr><td>R(1, 3)</td><td>-0.020</td><td>[-0.134;0.134]</td></tr><tr><td>R(2, 3)</td><td>0.158</td><td>[0.041;0.381]</td></tr></tbody></table>
# 6 Conclusion

The first aim of this paper was to provide a general framework for modeling multiple binary time series which can be compatible with many practical applications. In particular, we developed a general stationary theory as well as several asymptotic results for parametric estimation in multivariate autoregressive probit models. Our setup is very flexible since it allows to consider various kind of exogenous covariates, a single path or a longitudinal analysis and some pseudo-likelihood methods which are easy to implement numerically. A solid theoretical background is also provided and makes a bridge between some important results in probability theory and the asymptotics of pseudo-likelihood estimators for panel data.

Let us mention that our inference results could also be extended to other parametric models such as the multiple logistic model used for instance in Guanche et al. (2014) or Sebastian-Gonzalez et al. (2010). However, the use of a multivariate Gaussian distribution for the noise component $\varepsilon_t$ is a natural choice and allows to model many possible interactions between the binary coordinates. A more challenging extension concerns semi-parametric inference of the autoregressive parameters, when the probability distribution of the noise term $\varepsilon$ is not specified. In the case of binary choice models, corresponding to $k = 1$ in our context, Manski (1975) proposed the so-called maximum score estimator which was studied for a dynamic model by De Jong and Woutersen (2011). A multivariate version of such estimator could be interesting.

Another theoretical limitation of our study is the stationarity assumption imposed to the exogenous variables considered in our model. Indeed, many time series found for instance in macroeconomics or climate analysis cannot be considered as stationary. However, in a non-stationary framework, the asymptotic results could strongly depend on the dynamic of the covariates. See for instance Z. M. Debaly et al. (2023) for explosive covariates in count time series and also Park and Phillips (2000) for unit-root covariates in binary choice model.

Another way of improvement would be to extend the pairwise approach to composite likelihood obtained from sub-vectors of higher dimension such as triplets. We could then get more efficient estimators but a compromise between efficiency and numerical accuracy has to be done.

# 7 Supplementary material

## 7.1 Proof of Theorem 2.1

We first remind Theorem 3 from Truquet (2023).

**Theorem 7.1.** Let $(\zeta_t)_{t \in \mathbb{Z}}$ be a stationary and ergodic process valued in $\mathbb{R}^n$, and $E$ be a finite or countable set. For all $z \in \mathbb{R}^n$, we consider an application $F_z : E \to E$ such that the mapping

$$E \times \mathbb{R}^n \to E \\ (x, z) \mapsto F_z(x)$$

is measurable. For all $t, s$ in $\mathbb{Z}$ with $s \le t$, we denote the random mapping from $E$ to $E$

$$F_s^t = F_{\zeta_t} \circ \dots \circ F_{\zeta_s},$$

and we assume there exists $m \in \mathbb{N}^*$ such that

$$\mathbb{P}(\text{Card}(F_1^m(E)) = 1) > 0.$$

Then there exists a unique stationary and ergodic process $(Y_t)_{t \in \mathbb{Z}}$ such that

$$\forall t \in \mathbb{Z}, Y_t = F_{\zeta_t}(Y_{t-1}).$$

Furthermore, $(Y_t)_{t \in \mathbb{Z}}$ admits a Bernoulli shift representation

$$\forall t \in \mathbb{Z}, Y_t = G(\zeta_t, \zeta_{t-1}, \dots)$$

where $G : (\mathbb{R}^n)^N \to E$ is a measurable mapping.
1. We set $E = (\{0, 1\}^k)^p$ and we define the random mapping

$$F_{\zeta_t}(y_1, \dots, y_p) = (\mathbb{1}_{g(y_1, \dots, y_{p-1}, \zeta_t) > 0}, y_1, \dots, y_{p-1}),$$

where for $y \in \mathbb{R}^k$, $\mathbb{1}_{y>0}$ denotes the vector of $\mathbb{R}^k$ with binary coordinates $\mathbb{1}_{y_i>0}$, $1 \le i \le k$.

▷ We first prove that

$$\mathbb{P}(\text{Card}(F_1^p(E)) = 1) > 0.$$

To this end, note that

$$\bigcap_{i=1}^{k} C_{i}^{j_{i}} \subset \{\text{Card}(F_{1}^{p}(E)) = 1\}. \qquad (19)$$

Indeed, on the event $\cap_{i=1}^k C_i^{j_i}$, we have for any $(y_1, \dots, y_p) \in E$

$$F_{\zeta_1}(y_1, \dots, y_p) = (j, y_1, \dots, y_{p-1}),$$

where $j = (j_1, \dots, j_k)$. We then have

$$\begin{aligned} F_1^2(y_1, \dots, y_p) &= F_{\zeta_2}(j, y_1, \dots, y_{p-1}) \\ &= (\mathbb{1}_{g(j, y_1, \dots, y_{p-1}, \zeta_2) > 0}, j, y_1, \dots, y_{p-2}) \\ &= (j, j, y_1, \dots, y_{p-2}). \end{aligned}$$

Iterating this latter operation, we obtain $F_1^p(y_1, \dots, y_p) = (j, \dots, j)$, thus on the event $\cap_{i=1}^k C_i^{j_i}$, we have $F_1^p(E) = \{j\}$ and (19) follows.

By Assumption A2, we then have $\mathbb{P}(\text{Card}(F_1^p(E)) = 1) > 0$, and Theorem 7.1 gives the existence of an ergodic process $(Z_t)_{t \in \mathbb{Z}}$ such that

$$\forall t \in \mathbb{Z}, Z_t = F_{\zeta_t}(Z_{t-1}).$$

▷ We then define for all $t \in \mathbb{Z}$

$$Y_t = Z_t(1).$$

$(Y_t)_{t \in \mathbb{Z}}$ is an ergodic process and for all $t \in \mathbb{Z}$

$$Y_t = \mathbb{1}_{g(Z_{t-1}(1), \dots, Z_{t-1}(p), \zeta_t) > 0} = \mathbb{1}_{g(Y_{t-1}, Z_{t-1}(2), \dots, Z_{t-1}(p), \zeta_t) > 0}.$$

Yet, one can show by induction that for $t \in \mathbb{Z}$ and $l \in \{2, \dots, p\}$, $Z_t(l) = Z_{t-l+1}(1)$, thus

$$\forall t \in \mathbb{Z}, Y_t = \mathbb{1}_{g(Y'_{t-1}, \dots, Y'_{t-p}, \zeta_t) > 0}.$$

Furthermore, Theorem 7.1 specifies the existence of a measurable mapping

$$G : (\mathbb{R}^n)^N \to E$$

such that for all $t \in \mathbb{Z}$, $Z_t = G(\zeta_t, \zeta_{t-1}, \dots)$. Taking $H$ the composition $\pi_1 \circ G$, where $\pi_1$ denotes the projection on the first coordinate, we then obtain the Bernoulli shift representation (3).

▷ We finally prove the uniqueness of process $(Y_t)_{t \in \mathbb{Z}}$. Let $(Y'_t)_{t \in \mathbb{Z}}$ be a stationary and ergodic process satisfying

$$\forall t \in \mathbb{Z}, Y'_t = \mathbb{1}_{g(Y'_{t-1}, \dots, Y'_{t-p}, \zeta_t) > 0}$$

and denote $W_t = (Y'_t, \dots, Y'_{t-p+1})$. The process $(W_t)_{t \in \mathbb{Z}}$ is ergodic, and for all $t \in \mathbb{Z}$

$$\begin{aligned} F_{\zeta_t}(W_{t-1}) &= (\mathbb{1}_{g(Y'_{t-1}, \dots, Y'_{t-p}, \zeta_t) > 0}, Y'_{t-1}, \dots, Y'_{t-p+1}) \\ &= (Y'_{t-1}, \dots, Y'_{t-p+1}) \\ &= W_t. \end{aligned}$$
Since $(Z_t)_{t \in \mathbb{Z}}$ is the only process satisfying this latter recursion, we have

$$ \forall t \in \mathbb{Z}, W_t = Z_t, $$

and in particular

$$ \forall t \in \mathbb{Z}, Y_t = Y_t', $$

which proves the uniqueness of process $(Y_t)_{t \in \mathbb{Z}}$.

The process $(Y_t, \zeta_t)_{t \in \mathbb{Z}}$ is then ergodic, since we can write for all $t \in \mathbb{Z}$

$$ (Y_t, \zeta_t) = (H(\zeta_t, \zeta_{t-1}, \dots), \zeta_t). $$

2. We already established in the first point that $\cap_{i=1}^k C_i^{j_i} \subset \{\text{Card}(F_1^p(E) = 1)\}$, thus

$$ \mathbb{E} \left( \mathbb{1}_{\cap_{i=1}^k C_i^{j_i}} \mid \mathcal{F}_0 \right) \le \mathbb{E} \left( \mathbb{1}_{\text{Card}(F_1^p(E))=1} \mid \mathcal{F}_0 \right) $$

i.e.

$$ \mathbb{P} \left( \cap_{i=1}^{k} C_i^{j_i} \mid \mathcal{F}_0 \right) \le \mathbb{P} \left( \text{Card} (F_1^p(E)) = 1 \mid \mathcal{F}_0 \right), $$

hence

$$ \mathbb{P} (\text{Card} (F_1^p(E)) = 1 \mid \mathcal{F}_0) > 0 \text{ a.s.} $$

For simplicity of notations, we set $A_t = \{\text{Card}(F_{t-p+1}^p(E)) = 1\}$. We have

$$ \mathbb{1}_{A_t} = 1 \iff \exists y_0 \in E, \forall y \in E, F_{\zeta_t} \circ \dots \circ F_{\zeta_{t-p+1}}(y) = y_0. $$

For any $y \in E$, the application

$$ G_y : (\mathbb{R}^n)^p \to E \\ (e_1, \dots, e_p) \mapsto F_{e_p} \circ \dots \circ F_{e_1}(y) $$

is measurable, we thus have

$$ \mathbb{1}_{A_t} = \sum_{y_0 \in E} \prod_{y \in E} \mathbb{1}_{G_y^{-1}(\{y_0\})} (\zeta_t, \dots, \zeta_{t-p+1}) \\ = H(\zeta_t^{-}) $$

for some measurable function $H : (\mathbb{R}^n)^N \to \{0, 1\}$, where $\tau$ denotes the shift operator on the sequences.

For $t \in \mathbb{Z}$ fixed, we have by Birkhoff's ergodic theorem

$$ \frac{1}{S} \sum_{s=0}^{S-1} \mathbb{1}_{A_{t-s}} = \frac{1}{S} \sum_{s=0}^{S-1} H(\tau^s \zeta_t^{-}) \xrightarrow{S \to +\infty} \mathbb{P}(A_t \mid \mathcal{I}_t) \quad \text{a.s.} \quad (20) $$

where $\mathcal{I}_t = \{(\zeta_t^{-})^{-1}(I) : I \in \mathcal{B}(\mathbb{R}^s)^{\otimes N}, \tau^{-1}(I) = I\}$ and $\zeta_t^{-} = (\zeta_{t-j})_{j \ge 0}$.

Yet, by Lemma 7.2 below, we have $\mathcal{I}_t \subset \cap_{s \ge 0} \mathcal{F}_{t-s}$. In particular, $\mathcal{I}_t \subset \mathcal{F}_{t-p}$ and for all $J \in \mathcal{I}_t$

$$ \mathbb{P}(A_t \cap J) = \mathbb{E}[\mathbb{P}(A_t \mid \mathcal{I}_t) \mathbb{1}_J] = \mathbb{E}[\mathbb{P}(A_t \mid \mathcal{F}_{t-p}) \mathbb{1}_J]. $$

The choice $J = \{\mathbb{P}(A_t \mid \mathcal{I}_t) = 0\}$ then gives

$$ \mathbb{E}[\mathbb{P}(A_t \mid \mathcal{I}_t) \mathbb{1}_J] = 0, $$

and so $\mathbb{P}(J) = 0$ since $\mathbb{P}(A_t \mid \mathcal{F}_{t-p}) > 0$ almost surely.
Indeed, for all $K \in \mathcal{F}_{t-p}$, $K = (\zeta_{t-p}^{-})^{-1}(L)$, we have

$$
\begin{align*}
\mathbb{E}(\mathbb{1}_{A_t} \mathbb{1}_K) &= \mathbb{E}[\mathbb{P}(A_t | \mathcal{F}_{t-p}) \mathbb{1}_K] \\
&= \mathbb{E}[H(\zeta_t^{-}) \mathbb{1}_L(\zeta_{t-p}^{-})] \\
&= \mathbb{E}[H(\zeta_p^{-}) \mathbb{1}_L(\zeta_0^{-})] \quad \text{by stationarity} \\
&= \mathbb{E}[\mathbb{P}(A_p | \mathcal{F}_0) \mathbb{1}_L(\zeta_0^{-})].
\end{align*}
$$

The choice $K = \{\mathbb{P}(A_t | \mathcal{F}_{t-p}) = 0\}$ leads to

$$
\mathbb{E} [ \mathbb{P} ( A_p | \mathcal{F}_0 ) \mathbb{1}_L ( \zeta_0^- ) ] = 0,
$$

and since we have already proved that $\mathbb{P}(A_p | \mathcal{F}_0) > 0$ almost surely, we deduce

$$
\mathbb{P}(\zeta_0^- \in L) = \mathbb{P}(\zeta_{t-p}^- \in L) = \mathbb{P}(K) = 0,
$$

so $\mathbb{P}(A_t | \mathcal{F}_{t-p}) > 0$ almost surely.

It follows that $\mathbb{P}(A_t | \mathcal{I}_t) > 0$ almost surely, and by equation (20)

$$
\mathbb{P} \left( \limsup_{s \to +\infty} A_{t-s} \right) = 1.
$$

Moreover, if $\mathbb{1}_{A_{t-s}}(\omega) = 1$ for a certain value of $\omega$, we have for $T \ge p + s$

$$
\begin{align*}
F_{t-T}^{t}(y)(\omega) &= F_{t-s+1}^{t} \left( F_{t-s-p+1}^{t-s} \left( F_{t-T}^{t-s-p}(y) \right) \right) (\omega) \\
&= F_{t-s-p+1}^{t} (y_0)(\omega)
\end{align*}
$$

where $y_0 \in E$ is arbitrary. Then $\lim_{n \to \infty} F_{t-n}^t(y)$ exists a.s. and does not depend on $y$.

The end of the proof is similar to the one of Theorem 3 in Truquet (2023).

**Lemma 7.2.** For all $s \ge 0$, $\mathcal{I}_t \subset \mathcal{F}_{t-s}$.

*Proof.* We actually show by induction on $j$ that

$$
\forall j \in \mathbb{N}, \forall I \in \mathcal{B}(\mathbb{R}^s)^{\otimes \mathbb{N}}, \tau^{-1}(I) = I \implies (\zeta_t^{-})^{-1}(I) = (\zeta_{t-j}^{-})^{-1}(I).
$$

▷ For $j=0$, there is nothing to prove.

▷ Assume the property true for a fixed $j \in \mathbb{N}$. We thus have for all $I \in \mathcal{B}(\mathbb{R}^s)^{\otimes \mathbb{N}}$ such that $\tau^{-1}(I) = I$

$$
\begin{align*}
(\zeta_{t-j}^{-})^{-1}(I) &= \{\omega : (\zeta_{t-j}(\omega), \zeta_{t-j-1}(\omega), \dots) \in I\} \\
&= \{\omega : (\zeta_{t-j}(\omega), \zeta_{t-j-1}(\omega), \dots) \in \tau^{-1}I\} \\
&= \{\omega : (\zeta_{t-j-1}(\omega), \zeta_{t-j-2}(\omega), \dots) \in I\} \\
&= (\zeta_{t-j-1}^{-})^{-1}(I),
\end{align*}
$$

which concludes the induction.

<figure><img src="image_3.png" /></figure>
## 7.2 Proof of Corollary 2.2

Note first that $\sigma(X_s, \varepsilon_s | s \le t - 1) = \sigma(\zeta_{t-1}^{-}, X_{t-1})$.

The goal is to prove that assumption **A3** is satisfied. To this end, we set $j_1 = \dots = j_k = 1$ and for $t \in \{1, \dots, p\}$

$$B_t = \bigcap_{i=1}^{k} \left\{ \varepsilon_{i,t} > - \min_{y_1, \dots, y_p} h_i(y_1, \dots, y_p, X_{t-1}) \right\},$$

so that we have $\cap_{i=1}^k C_i^{j_i} = \cap_{t=1}^p B_t$.

By Assumption **A4**, we have

$$\mathbb{P}(B_t | \sigma(\zeta_{t-1}^{-}, X_{t-1})) > 0 \text{ a.s.}$$

Indeed, for any $C \in \sigma(\zeta_{t-1}^{-}, X_{t-1})$, $C = (\zeta_{t-1}^{-}, X_{t-1})^{-1}(B)$, we have

$$\begin{aligned} \mathbb{E}(\mathbb{1}_{B_t}\mathbb{1}_C) &= \mathbb{E}\left[\prod_{i=1}^{k} \mathbb{1}\left\{\varepsilon_{i,t} > - \min_{y_1, \dots, y_p} h_i(y_1, \dots, y_p, X_{t-1})\right\}\mathbb{1}_C\right] \\ &= \mathbb{E}(\mathbb{1}_K(X_{t-1}, \varepsilon_t)\mathbb{1}_C) \end{aligned}$$

where $K = \{(x, e) \in \mathbb{R}^d \times \mathbb{R}^k : \forall i \in \{1, \dots, k\}, x_i > - \min_{y_1, \dots, y_p} h_i(y_1, \dots, y_p, e)\}$. Thus

$$\begin{aligned} \mathbb{E}(\mathbb{1}_{B_t}\mathbb{1}_C) &= \int \mathbb{1}_K(X_{t-1}(\omega), \varepsilon_t(\omega)) \mathbb{1}_B(\zeta_{t-1}^{-}(\omega), X_{t-1}(\omega)) \mathbb{P}(d\omega) \\ &= \int \mathbb{1}_K(x_{t-1}, e) \mathbb{1}_B(s^{-}, x_{t-1}) \mathbb{P}_{(\varepsilon_t, \zeta_{t-1}^{-}, X_{t-1})}(de, ds^{-}, dx_{t-1}) \\ &= \iint \mathbb{1}_K(x_{t-1}, e) \mathbb{P}_{\varepsilon_t|(\zeta_{t-1}^{-}, X_{t-1})}(de | s^{-}, x_{t-1}) \mathbb{1}_B(s^{-}, x_{t-1}) \mathbb{P}_{(\zeta_{t-1}^{-}, X_{t-1})}(ds^{-}, dx_{t-1}) \end{aligned}$$

where $\mathbb{P}_{\varepsilon_t|(\zeta_{t-1}^{-}, X_{t-1})}$ is a conditional distribution of $\varepsilon_t$ given $(\zeta_{t-1}^{-}, X_{t-1})$

$$= \int \mathbb{P}_{\varepsilon_t|(\zeta_{t-1}^{-}, X_{t-1})}(K_{x_{t-1}} | s^{-}, x_{t-1}) \mathbb{1}_B(s^{-}, x_{t-1}) \mathbb{P}_{(\zeta_{t-1}^{-}, X_{t-1})}(ds^{-}, dx_{t-1})$$

with

$$K_{x_{t-1}} = \left\{ e \in \mathbb{R}^k : \forall i \in \{1, \dots, k\}, e_i > - \min_{y_1, \dots, y_p} h_i(y_1, \dots, y_p, x_{t-1}) \right\} = \cap_{i=1}^k \varphi_{i,x_{t-1}}^{-1}((0, +\infty))$$

and $\varphi_{i,x_{t-1}}(e) = e_i + \min_{y_1, \dots, y_p} h_i(y_1, \dots, y_p, x_{t-1})$ is trivially continuous. In particular, $K_{x_{t-1}}$ is a non-empty open set, and by assumption **A4** we have $\mathbb{P}_{(\zeta_{t-1}^{-}, X_{t-1})}$-almost surely

$$\mathbb{P}_{\varepsilon_t|(\zeta_{t-1}^{-}, X_{t-1})}(K_{x_{t-1}} | s^{-}, x_{t-1}) > 0.$$

If we set $C = \{\omega : \mathbb{P}(B_t | \sigma(\zeta_{t-1}^{-}, X_{t-1})) (\omega) = 0\}$, and suppose by absurd that $\mathbb{P}(C) > 0$, we then get the contradiction

$$0 = \mathbb{E}(\mathbb{1}_{B_t}\mathbb{1}_C) = \int \mathbb{P}_{\varepsilon_t|(\zeta_{t-1}^{-}, X_{t-1})}(K_{x_{t-1}} | s^{-}, x_{t-1}) \mathbb{1}_B(s^{-}, x_{t-1}) \mathbb{P}_{(\zeta_{t-1}^{-}, X_{t-1})}(ds^{-}, dx_{t-1}) > 0.$$

Thus $\mathbb{P}(C) = 0$, i.e.

$$\mathbb{P}(B_t | \sigma(\zeta_{t-1}^{-}, X_{t-1})) > 0 \text{ a.s.} \quad (21)$$

We then set for $j \in \{1, \dots, p\}$

$$U_j = \mathbb{P}\left(\cap_{t=1}^{j} B_t \mid \mathcal{F}_0\right),$$

and we show by induction on $j$ that

$$\forall j \in \{1, \dots, p\}, U_j > 0 \text{ a.s.}$$
▷ We already have $U_1 > 0$ almost surely by equation (21).

▷ Assume that $U_j > 0$ almost surely for a fixed $j \in \{1, \dots, p\}$, and take $C \in \mathcal{F}_0$, we have

$$ \begin{aligned} \mathbb{E} (\mathbb{1}_{\cap_{t=1}^{j+1} B_j} \mathbb{1}_C) &= \mathbb{E} (\mathbb{1}_{B_{j+1}} \mathbb{1}_{\cap_{j=1}^t B_j} \mathbb{1}_C) \\ &= \mathbb{E} (\mathbb{P}(B_{j+1} | \mathcal{F}_j) \mathbb{1}_{\cap_{j=1}^t B_j} \mathbb{1}_C), \end{aligned} $$

thus $U_{j+1} = \mathbb{E} (\mathbb{P}(B_{j+1} | \mathcal{F}_j) \mathbb{1}_{\cap_{j=1}^t B_j} | \mathcal{F}_0)$. However, $\mathbb{P}(B_{j+1} | \mathcal{F}_j) > 0$ almost surely by equation (21), and $\mathbb{E} (\mathbb{1}_{\cap_{t=1}^j B_t} | \mathcal{F}_0) > 0$ almost surely by assumption. Thus we have $U_{j+1} > 0$ almost surely, which concludes the induction.

We conclude that

$$ U_p = \mathbb{P}(\cap_{t=1}^p B_t | \mathcal{F}_0) = \mathbb{P}(\cap_{i=1}^k C_i^{j_i} | \mathcal{F}_0) > 0 \quad \text{a.s.} $$

and assumption **A3** is verified.

# 7.3 Proof of Theorem 3.2

## 7.3.1 Proof for estimator $\hat{\theta}_{\mathcal{P}\mathcal{L}}$

The consistency, the stochastic expansion, and the asymptotic normality of the estimator $\hat{\theta}$ is obtained by using Theorem 7.9 and Theorem 7.10. We thus have to check the assumptions **H1** to **H9** mentioned in these latter theorems for the process $(Z_t)_{t \in \mathbb{Z}}$ defined by

$$ Z_t = (Y_t, \dots, Y_{t-p}, X_{t-1}), $$

and the family of applications $(m_\theta)_{\theta \in \Theta}$ defined by

$$ m_{\theta}(Z_t) = \log \left( \int_{\mathbb{R}^k} \prod_{i=1}^k \mathbb{1}_{I_{Y_i, t}} (\lambda_{i,t} + e_i) \varphi_R(e) \mathrm{d}e \right). $$

▷ **H1** is satisfied since $(Y_t)_{t \in \mathbb{Z}}$ and $(X_t)_{t \in \mathbb{Z}}$ are both ergodic by assumption **B1** and corollary 2.2.

▷ **H2** is satisfied by assumption **B4**.

▷ In order to check assumption **H3**, we first prove the continuity of the application $\theta \mapsto m_\theta(Z_0)$ for any fixed value of $Z_0$. Here we have

$$ m_{\theta}(Z_0) = \log(f_{\theta}(Z_0)) $$

where

$$ f_{\theta}(Z_0) = \int_{\mathbb{R}^k} \prod_{i=1}^k \mathbb{1}_{I_{Y_{i,0}} - \lambda_{i,0}(\gamma)}(x_i) \varphi_R(x) \mathrm{d}x. $$

Let us fix $\theta_0 = (\gamma_0, R_0) \in \Theta$, for any $i \in \{1, \dots, k\}$, we denote

$$ H_i = \{x = (x_1, \dots, x_k) \in \mathbb{R}^k : x_i = -\lambda_{i,0}(\gamma_0)\} $$

and $H = \cup_{i=1}^k H_i$. $H$ is negligible for the Lebesgue's measure, since it is a finite union of hyperplanes, and it is easy to verify that

$$ \prod_{i=1}^{k} \mathbb{1}_{I_{Y_{i,0}} - \lambda_{i,0}(\gamma)}(x_i) \varphi_R(x) \xrightarrow{\theta \to \theta_0} \prod_{i=1}^{k} \mathbb{1}_{I_{Y_{i,0}} - \lambda_{i,0}(\gamma_0)}(x_i) \varphi_{R_0}(x) $$

for any $x \notin H$.
Furthermore, according to Lemma 7.3 below, it is possible to find positive constants $m_{\det}$ and $c_2$ such that for any $\theta \in \Theta$ and any $x \in \mathbb{R}^k$

$$ \prod_{i=1}^{k} 1_{I_{Y_{i,0}} - \lambda_{i,0}(\gamma)}(x_i) \varphi_R(x) \le \frac{1}{\sqrt{(2\pi)^k m_{\det}}} \exp \left( -\frac{c_2}{2} \|x\|_2^2 \right), $$

which is integrable on $\mathbb{R}^k$. Thus, by continuity under the integral sign, the application $\theta \mapsto f_{\theta}(Z_0)$ is continuous, and it follows that $\theta \mapsto m_{\theta}(Z_0)$ is also continuous.

In order prove that $\mathbb{E}(\sup_{\theta \in \Theta} |m_{\theta}(Z_0)|) < +\infty$, we use once again Lemma 7.3 to obtain positive constants $M_{\det}$ and $c_1$ such that for any $\theta \in \Theta$ and $x \in \mathbb{R}^k$

$$ \begin{aligned} f_{\theta}(Z_0) &\ge \frac{1}{\sqrt{(2\pi)^k M_{\det}}} \int_{\mathbb{R}^k} \prod_{i=1}^k 1_{I_{Y_{i,0}} - \lambda_{i,0}(\gamma)}(x_i) \exp \left( -\frac{c_1}{2} \|x\|_2^2 \right) dx \\ &\ge \frac{1}{\sqrt{c_1 M_{\det}}} \prod_{i=1}^k \min \{ \Phi(-\sqrt{c_1} \lambda_{i,0}(\gamma)), 1 - \Phi(-\sqrt{c_1} \lambda_{i,0}(\gamma)) \} \end{aligned} $$

$$ \begin{aligned} \text{with the substitution } &u = \sqrt{c_1} x \\ &\ge \frac{1}{\sqrt{c_1 M_{\det}}} \exp \left( -dc_1 \|\lambda_0(\gamma)\|_2^2 \right) \end{aligned} $$

where $d$ is the positive constant mentioned in Lemma 7.8. For any $\theta \in \Theta$, we thus have

$$ \begin{aligned} |m_{\theta}(Z_0)| &= -\log(f_{\theta}(Z_0)) \\ &\le \frac{1}{2} \log(c_1 M_{\det}) + dc_1 \|\lambda_0(\gamma)\|_2^2 \\ &\le \frac{1}{2} \log(c_1 M_{\det}) + dc_1 (M_A + M_B \|X_{-1}\|_2)^2 \end{aligned} $$

for some positive constants $M_A$ and $M_B$ obtained by compactness of $\Theta$. Since it is assumed that $X_{-1} \in L^2$, this latter quantity is integrable and we have

$$ \mathbb{E} \left( \sup_{\theta \in \Theta} |m_{\theta}(Z_0)| \right) < +\infty. $$

**Lemma 7.3.** There exist constants $M_{\det}$, $m_{\det}$, $c_1$, $c_2 > 0$ such that for all $R \in \mathbb{R}$ and all $x \in \mathbb{R}^k$

$$ \frac{1}{\sqrt{(2\pi)^k M_{\det}}} \times \exp \left( -\frac{c_1}{2} \|x\|_2^2 \right) \le \varphi_R(x) \le \frac{1}{\sqrt{(2\pi)^k m_{\det}}} \times \exp \left( -\frac{c_2}{2} \|x\|_2^2 \right). $$

*Proof of Lemma 7.3.* By continuity of the determinant, and by compactness of $\mathbb{R}$, there exist $m_{\det}$, $M_{\det} > 0$ such that

$$ \forall R \in \mathbb{R}, m_{\det} \le \det(R) \le M_{\det}. $$

Furthermore, take any $R \in \mathbb{R}$, we have for all $x \in \mathbb{R}^k$

$$ \begin{aligned} \varphi_R(x) &= \frac{1}{\sqrt{(2\pi)^k \det(R)}} \exp \left( -\frac{1}{2} \langle R^{-1}x, x \rangle \right) \\ &= \frac{1}{\sqrt{(2\pi)^k \det(R)}} \exp \left( -\frac{1}{2} \|x\|_{R^{-1}}^2 \right) \end{aligned} \qquad (22) $$

where $\| \cdot \|_{R^{-1}}$ denotes the norm associated with the matrix $R^{-1}$. By equivalence of the norms, there exist $c(R), C(R) > 0$ such that

$$ \forall x \in \mathbb{R}^k, c(R) \times \|x\|_2 \le \|x\|_{R^{-1}} \le C(R) \times \|x\|_2. $$
Without loss of generality, we can suppose that

$$c(R) = \sup \left\{ c > 0 : \forall x \in \mathbb{R}^k, c \|x\|_2 \le \|x\|_{R^{-1}} \right\}$$

and

$$C(R) = \inf \left\{ C > 0 : \forall x \in \mathbb{R}^k, C \|x\|_2 \ge \|x\|_{R^{-1}} \right\}.$$

Let us show that there exists $c > 0$ such that

$$\forall R \in \mathcal{R}, c \le c(R).$$

By absurd, assume that for all integer $n \ge 1$, there exists $R_n \in \mathcal{R}$ such that

$$c(R_n) < \frac{1}{n}.$$

Since $\mathcal{R}$ is compact, it is possible to extract a subsequence $(R_{\varphi(n)})_{n \ge 1}$ that converges toward $\tilde{R} \in \mathcal{R}$.
Yet, for all $n \ge 1$, we have $c(R_{\varphi(n)}) < \frac{1}{\varphi(n)}$, thus there exists $x_{\varphi(n)} \in \mathcal{S}^1$ such that

$$\|x_{\varphi(n)}\|_{R^{-1}} < \frac{1}{\varphi(n)}.$$

By compactness once again, we can extract a subsequence $(x_{\varphi \circ \psi(n)})_{n \ge 1}$ that converges toward $\tilde{x} \in \mathcal{S}^1$.
Thus, for all $n \ge 1$

$$\|x_{\varphi \circ \psi(n)}\|_{R^{-1}} < \frac{1}{\varphi \circ \psi(n)} \xrightarrow{n \to +\infty} 0.$$

However, by continuity of the inverse, we also have $\|x_{\varphi \circ \psi(n)}\|_{R^{-1}} \xrightarrow{n \to +\infty} \|\tilde{x}\|_{\tilde{R}^{-1}}$. Hence $\tilde{x} = 0$, which is a contradiction.

Therefore, there exists indeed a constant $c > 0$ such that

$$\forall R \in \mathcal{R}, c \le c(R).$$

One can show the same way that there exists a constant $C > 0$ such that

$$\forall R \in \mathcal{R}, C(R) \le C.$$

Combining with equality (22), we finally have

$$\frac{1}{\sqrt{(2\pi)^k M_{\det}}} \times \exp\left(-\frac{C^2}{2} \|x\|_2^2\right) \le \varphi_R(x) \le \frac{1}{\sqrt{(2\pi)^k m_{\det}}} \times \exp\left(-\frac{c^2}{2} \|x\|_2^2\right),$$

which is the desired result with $c_1 = C^2$ and $c_2 = c^2$. $\square$

▷ We now check **H4**. For $\theta \in \Theta$ fixed, setting for any $t \in \mathbb{Z}$, $Y_t^- = (Y_t, Y_{t-1}, \dots)$, we have

$$\mathbb{E}(m_{\theta}(Z_0)) - \mathbb{E}(m_{\theta_0}(Z_0)) = \mathbb{E}\left(\log\left(\frac{\psi_{\theta}(Y_0, Y_{-1}^{-}, X_{-1})}{\psi_{\theta_0}(Y_0, Y_{-1}^{-}, X_{-1})}\right)\right),$$

where by assumption **B2**, $\psi_{\theta}(y_0, Y_{-1}^{-}, X_{-1}) = \int_{\mathbb{R}^k} \prod_{i=1}^k 1_{I_{y_{i,0}}}(\lambda_{i,0}(\gamma) + x_i)\varphi_R(x)dx$ is the density of the conditional distribution $\mathcal{L}_{\theta}(Y_0 | Y_{-1}^{-}, X_{-1})$ with respect to the measure $\mu = \sum_{y \in \{0,1\}^k} \delta_y$ on $\{0,1\}^k$.

Using the fact that

$$\forall x > 0, \log(x) \le 2(\sqrt{x} - 1),$$
with equality if and only if $x = 1$, we thus have

$$
\begin{align*}
\mathbb{E}(m_{\theta}(Z_0)) - \mathbb{E}(m_{\theta_0}(Z_0)) &= \iint \log \left( \frac{\psi_{\theta}(y_0, y_{-1}^{-}, x_{-1})}{\psi_{\theta_0}(y_0, y_{-1}^{-}, x_{-1})} \right) \psi_{\theta_0}(y_0, y_{-1}^{-}, x_{-1}) \mu(dy_0) \mathbb{P}_{(Y_{-1}^{-}, X_{-1})}(dy_{-1}^{-}, dx_{-1}) \\
&\le \iint 2\sqrt{\psi_{\theta}(y_0, y_{-1}^{-}, x_{-1})} \times \sqrt{\psi_{\theta_0}(y_0, y_{-1}^{-}, x_{-1})} \mu(dy_0) - 2\mathbb{P}_{(Y_{-1}^{-}, X_{-1})}(dy_{-1}^{-}, dx_{-1}) \\
&= -\iint \left( \sqrt{\psi_{\theta}(y_0, y_{-1}^{-}, x_{-1})} - \sqrt{\psi_{\theta_0}(y_0, y_{-1}^{-}, x_{-1})} \right)^2 \mu(dy_0) \mathbb{P}_{(Y_{-1}^{-}, X_{-1})}(dy_{-1}^{-}, dx_{-1}) \\
&\le 0,
\end{align*}
$$

hence

$$
\forall \theta \in \Theta, \mathbb{E}(m_{\theta}(Z_0)) \le \mathbb{E}(m_{\theta_0}(Z_0)).
$$

Furthermore, the equality case implies that $\mathbb{P}_{(Y_{-1}^{-}, X_{-1})}$-almost surely

$$
\int \left[ \log \left( \frac{\psi_{\theta}(y_0, y_{-1}^{-}, x_{-1})}{\psi_{\theta_0}(y_0, y_{-1}^{-}, x_{-1})} \right) - 2 \left( \sqrt{\frac{\psi_{\theta}(y_0, y_{-1}^{-}, x_{-1})}{\psi_{\theta_0}(y_0, y_{-1}^{-}, x_{-1})}} - 1 \right) \right] \psi_{\theta_0}(y_0, y_{-1}^{-}, x_{-1}) \mu(dy_0) = 0,
$$

thus, $\mathbb{P}_{(Y_{-1}^{-}, X_{-1})}$-almost surely, we have for any $y_0 \in \{0, 1\}^k$

$$
\left[ \log \left( \frac{\psi_{\theta}(y_0, y_{-1}^{-}, x_{-1})}{\psi_{\theta_0}(y_0, y_{-1}^{-}, x_{-1})} \right) - 2 \left( \sqrt{\frac{\psi_{\theta}(y_0, y_{-1}^{-}, x_{-1})}{\psi_{\theta_0}(y_0, y_{-1}^{-}, x_{-1})}} - 1 \right) \right] \psi_{\theta_0}(y_0, y_{-1}^{-}, x_{-1}) = 0.
$$

Since $\psi_{\theta_0}(y_0, y_{-1}^{-}, x_{-1})$ is necessarily positive, whatever the values of $y_0, y_{-1}^{-}$ and $x_{-1}$, we have $\mathbb{P}_{(Y_{-1}^{-}, X_{-1})}$-almost surely

$$
\forall y_0 \in \{0, 1\}^k, \log \left( \frac{\psi_{\theta}(y_0, y_{-1}^{-}, x_{-1})}{\psi_{\theta_0}(y_0, y_{-1}^{-}, x_{-1})} \right) - 2 \left( \sqrt{\frac{\psi_{\theta}(y_0, y_{-1}^{-}, x_{-1})}{\psi_{\theta_0}(y_0, y_{-1}^{-}, x_{-1})}} - 1 \right) = 0,
$$

which leads to

$$
\forall y_0 \in \{0, 1\}^k, \quad \psi_{\theta}(y_0, y_{-1}^{-}, x_{-1}) = \psi_{\theta_0}(y_0, y_{-1}^{-}, x_{-1})
$$

$\mathbb{P}_{(Y_{-1}^{-}, X_{-1})}$-almost surely. Equivalently, this latter equation means that the conditional distributions
$\mathcal{L}_{\theta}(Y_0, | Y_{-1}^{-}, X_{-1})$ and $\mathcal{L}_{\theta_0}(Y_0, | Y_{-1}^{-}, X_{-1})$ are almost surely identical.

In particular, we have almost surely

$$
\forall i \in \{1, \dots, k\}, \quad \forall y \in \{0, 1\}, \quad \mathbb{P}_{\theta} (Y_{i,0} = y | Y_{-1}^{-}, X_{-1}) = \mathbb{P}_{\theta_0} (Y_{i,0} = y | Y_{-1}^{-}, X_{-1}),
$$

thus we have almost surely

$$
\forall i \in \{1, \dots, k\}, \quad \mathbb{P}_{\theta} (Y_{i,0} = 1 | Y_{-1}^{-}, X_{-1}) = \mathbb{P}_{\theta_0} (Y_{i,0} = 1 | Y_{-1}^{-}, X_{-1})
$$

$$
\text{i.e. } \forall i \in \{1, \dots, k\}, \quad \Phi(-\lambda_{i,0}(\gamma)) = \Phi(-\lambda_{i,0}(\gamma_0))
$$

$$
\text{i.e. } \forall i \in \{1, \dots, k\}, \quad \lambda_{i,0}(\gamma) = \lambda_{i,0}(\gamma_0).
$$

If we denote $A_1^{(0)}, \dots, A_p^{(0)}, B^{(0)}$ and $R^{(0)}$ the values of the matrices composing $\theta_0$, this latter equation
leads to

$$
\sum_{l=1}^{p} \underbrace{(A_l - A_l^{(0)})}_{M_l} Y_{-l} + \underbrace{(B - B^{(0)})}_{N} X_{-1} = 0 \quad \text{a.s.} \qquad (23)
$$

By assumption **B5**, we immediately have $N = 0$.

Furthermore, since for any $y_{-1}, \dots, y_{-p}$ in $\{0, 1\}^k$,

$$
\mathbb{P}(\{(Y_{-1}, \dots, Y_{-p}) = (y_{-1}, \dots, y_{-p})\}) > 0,
$$
we have for any $y_{-1}, \dots, y_{-p}$ in $\{0, 1\}^k$

$$M_1 \cdot y_{-1} = - \sum_{l=2}^{p} M_l \cdot y_{-l}.$$

This equality is true in the particular case where $y_{-2} = \dots = y_{-p} = 0$, and then we have

$$\forall y_{-1} \in \{0, 1\}^k, \quad M_1 \cdot y_{-1} = 0.$$

This latter equality is thus true for any vector in the canonical base of $\mathbb{R}^k$, and it follows that $M_1 = 0$. A simple induction then leads to

$$\forall l \in \{1, \dots, p\}, \quad M_l = 0.$$

**Note.** In the specific case where the processes $(X_t)_{t \in \mathbb{Z}}$ and $(\varepsilon_t)_{t \in \mathbb{Z}}$ are supposed independent, and assumption **B5** is replaced by assumption **B5'**, the previous proof can be modified from equality (23). In this case, let us suppose by absurd that $M_1 \neq 0$. Then, there exists $v \in \mathbb{R}^k \setminus \{0\}$ such that

$$1 = \mathbb{P}(v' \cdot Y_{-1} = G(X_{-1}^{-}, \varepsilon_{-2}^{-}))$$

where $G$ is a measurable application

$$= \iint \mathbb{1}_{\{G(x_{-1}^{-}, e_{-2}^{-})\}} (v' \cdot y_{-1}) \mathbb{P}_{Y_{-1}|(X_{-1}^{-}, \varepsilon_{-2}^{-})} (x_{-1}^{-}, e_{-2}^{-}) (dy_{-1}) \mathbb{P}_{(X_{-1}^{-}, \varepsilon_{-2}^{-})} (dx_{-1}^{-}, de_{-2}^{-}).$$

Thus, for $\mathbb{P}_{(X_{-1}^{-}, \varepsilon_{-2}^{-})}$-almost all $(x_{-1}^{-}, e_{-2}^{-})$ we have

$$\int \mathbb{1}_{\{G(x_{-1}^{-}, e_{-2}^{-})\}} (v' \cdot y_{-1}) \mathbb{P}_{Y_{-1}|(X_{-1}^{-}, \varepsilon_{-2}^{-})} (x_{-1}^{-}, e_{-2}^{-}) (dy_{-1}) = 1.$$

There exists in particular an element $(x_{-1}^{-}, e_{-2}^{-})$ such that for $\mathbb{P}_{Y_{-1}|(X_{-1}^{-}, \varepsilon_{-2}^{-})} (x_{-1}^{-}, e_{-2}^{-})$-almost all $y \in \{0, 1\}^k$,

$$v' \cdot y = G(x_{-1}^{-}, e_{-2}^{-}). \qquad (24)$$

However, the conditional distribution $\mathcal{L}_{\theta}(Y_{-1} | X_{-1}^{-}, \varepsilon_{-2}^{-})$ is fully supported by $\{0, 1\}^k$. Denoting $H(X_{-1}^{-}, \varepsilon_{-2}^{-}) = \sum_{l=1}^p A_l Y_{-l-1} + B X_{-2}$, we have indeed for any adequate measurable set $C$ and any $y_{-1} \in \{0, 1\}^k$

$$\begin{align*}
\mathbb{P}(Y_{-1} = y_{-1}, (X_{-1}^{-}, \varepsilon_{-2}^{-}) \in C) &= \mathbb{P} \left( \varepsilon_{-1} \in \prod_{i=1}^k I_{y_{i,-1}-H(X_{-1}^{-}, \varepsilon_{-2}^{-})}, (X_{-1}^{-}, \varepsilon_{-2}^{-}) \in C \right) \\
&= \iint_{\mathbb{R}^k} \prod_{i=1}^k \mathbb{1}_{I_{y_{i,-1}-H(x_{-1}^{-}, e_{-2}^{-})} (e_{-1})} \mathbb{P}_{\varepsilon_{-1}|(X_{-1}^{-}, \varepsilon_{-2}^{-})} (x_{-1}^{-}, e_{-2}^{-}) (de_{-1}) \\
&\qquad \times \mathbb{1}_C(x_{-1}^{-}, e_{-2}^{-}) \mathbb{P}_{(X_{-1}^{-}, \varepsilon_{-2}^{-})} (dx_{-1}^{-}, de_{-2}^{-}) \\
&= \iint_{\mathbb{R}^k} \prod_{i=1}^k \mathbb{1}_{I_{y_{i,-1}-H(x_{-1}^{-}, e_{-2}^{-})} (e_{i,-1})} \mathbb{P}_{\varepsilon_{-1}} (de_{-1}) \\
&\qquad \times \mathbb{1}_C(x_{-1}^{-}, e_{-2}^{-}) \mathbb{P}_{(X_{-1}^{-}, \varepsilon_{-2}^{-})} (dx_{-1}^{-}, de_{-2}^{-})
\end{align*}$$

since $\varepsilon_{-1}$ is independent from $(X_{-1}^{-}, \varepsilon_{-2}^{-})$ from the assumption of independence between $(X_t)_{t \in \mathbb{Z}}$ and $(\varepsilon_t)_{t \in \mathbb{Z}}$ and assumptions **B2**. Thus

$$\mathbb{P}(Y_{-1} = y_{-1} | (X_{-1}^{-}, \varepsilon_{-2}^{-}) = (x_{-1}^{-}, e_{-2}^{-})) = \int_{\mathbb{R}^k} \prod_{i=1}^k \mathbb{1}_{I_{y_{i,-1}-G(x_{-1}^{-}, e_{-2}^{-})} (e_{i,-1})} \mathbb{P}_{\varepsilon_{-1}} (de_{-1}) > 0.$$
Choosing $y = 0$ in (24), we obtain $G(x_{-1}^{-}, e_{-2}^{-}) = 0$, and it follows that for any $y \in \{0, 1\}^k$

$$v' \cdot y = 0,$$

which is impossible since $v \neq 0$. Thus, $M_1 = 0$, and an immediate induction leads to

$$\forall l \in \{1, \dots, p\}, \quad M_l = 0.$$

Equation (23) then becomes

$$N \cdot X_{-1} = 0,$$

and assumption **B5'** gives $N = 0$.

In the end, we have shown that the equality $\mathbb{E}(m_{\theta}(Z_0)) = \mathbb{E}(m_{\theta_0}(Z_0))$ implies $B = B^{(0)}$ and $A_l^{(0)} = A_l$ for any $l \in \{1, \dots, p\}$.

We still have to prove that this latter equality also implies that $R = R^{(0)}$. We already mentioned that the conditional distributions $\mathcal{L}_{\theta}(Y_0 | Y_{-1}^{-}, X_{-1})$ and $\mathcal{L}_{\theta_0}(Y_0 | Y_{-1}^{-}, X_{-1})$ are almost surely identical. Therefore, it is also the case for the distributions of the couples, i.e. for any $i < j$ in $\{1, \dots, k\}$

$$\mathcal{L}_{\theta}((Y_{i,0}, Y_{j,0}) | Y_{-1}^{-}, X_{-1}) \quad \text{and} \quad \mathcal{L}_{\theta_0}((Y_{i,0}, Y_{j,0}) | Y_{-1}^{-}, X_{-1})$$

are almost surely identical. Note that if $k = 1$, there is nothing to prove since the covariance matrix $R$ is actually a correlation matrix.

In particular, we have almost surely

$$\mathbb{P}_{\theta} (Y_{i,0} = 0, Y_{j,0} = 0 | Y_{-1}^{-}, X_{-1}) = \mathbb{P}_{\theta_0} (Y_{i,0} = 0, Y_{j,0} = 0 | Y_{-1}^{-}, X_{-1})$$

i.e.

$$\int_{-\infty}^{-\lambda_{i,0}(\gamma)} \Phi \left( \frac{-\lambda_{j,0}(\gamma) - R(i, j)x_i}{\sqrt{1 - R(i, j)^2}} \right) \varphi(x_i) dx_i = \int_{-\infty}^{-\lambda_{i,0}(\gamma_0)} \Phi \left( \frac{-\lambda_{j,0}(\gamma_0) - R^{(0)}(i, j)x_i}{\sqrt{1 - R^{(0)}(i, j)^2}} \right) \varphi(x_i) dx_i.$$

Since we have already proved that $\lambda_0(\gamma) = \lambda_0(\gamma_0)$, we thus have almost surely

$$\int_{-\infty}^{-\lambda_{i,0}(\gamma_0)} \Phi \left( \frac{-\lambda_{j,0}(\gamma_0) - R(i, j)x_i}{\sqrt{1 - R(i, j)^2}} \right) \varphi(x_i) dx_i = \int_{-\infty}^{-\lambda_{i,0}(\gamma_0)} \Phi \left( \frac{-\lambda_{j,0}(\gamma_0) - R^{(0)}(i, j)x_i}{\sqrt{1 - R^{(0)}(i, j)^2}} \right) \varphi(x_i) dx_i. \quad (25)$$

However, for any $\omega_i, \omega_j$ in $\mathbb{R}$, the mapping

$$f: \begin{array}{r} (-1, 1) \rightarrow \mathbb{R} \\ r \mapsto \int_{-\infty}^{\omega_i} \Phi \left( \frac{\omega_j - r x_i}{\sqrt{1 - r^2}} \right) \varphi(x_i) dx_i \end{array}$$

is an increasing function (see Lemma 7.4). We then deduce that $R(i, j) = R^{(0)}(i, j)$ for any $i < j$ in $\{1, \dots, k\}$, thus $R = R^{(0)}$.

As a consequence, we do have

$$\forall \theta \in \Theta, \mathbb{E}(m_{\theta}(Z_0)) \leq \mathbb{E}(m_{\theta_0}(Z_0)),$$

with equality if and only if $\theta = \theta_0$, so **H4** is satisfied.

**Lemma 7.4.** Take any $\omega_i, \omega_j$ in $\mathbb{R}$, the function

$$f: \begin{array}{r} (-1, 1) \rightarrow \mathbb{R} \\ r \mapsto \int_{-\infty}^{\omega_i} \Phi \left( \frac{\omega_j - r e_i}{\sqrt{1 - r^2}} \right) \varphi(e_i) de_i \end{array}$$

is strictly increasing.
Proof of Lemma 7.4. We first apply the theorem of differentiation of a function defined by an integral to $f$. To do so, we verify the following assertions.

▷ For any $r \in (-1, 1)$ and any $x_i \in (-\infty, \omega_i]$, we have

$$ \left| \Phi \left( \frac{\omega_j - rx_i}{\sqrt{1-r^2}} \right) \varphi(x_i) \right| \le \varphi(x_i), $$

thus the integrand function defining $f$ is integrable on $]-\infty, \omega_i]$.

▷ For any $x_i \in (-\infty, \omega_i]$, the application

$$ r \mapsto \Phi \left( \frac{\omega_j - rx_i}{\sqrt{1-r^2}} \right) \varphi(x_i) $$

is differentiable, and its derivative is equal to

$$ \frac{r\omega_j - x_i}{(1-r^2)^{3/2}} \times \varphi \left( \frac{\omega_j - rx_i}{\sqrt{1-r^2}} \right) \varphi(x_i). $$

Without loss of generality, we can take $r$ in a compact subset of $(-1, 1)$, and this latter function can be bounded by an integrable function.

The function $f$ is thus differentiable on $(-1, 1)$ and we have

$$ \forall r \in (-1, 1), \quad \dot{f}(r) = \int_{-\infty}^{\omega_i} \frac{r\omega_j - x_i}{(1-r^2)^{3/2}} \times \varphi\left(\frac{\omega_j - rx_i}{\sqrt{1-r^2}}\right) \varphi(x_i)dx_i. $$

Denoting $\alpha = \frac{\omega_j}{\sqrt{1-r^2}}$, $\beta = \frac{r}{\sqrt{1-r^2}}$, $\gamma = \frac{r\alpha}{1-r^2}$ and $\delta = \frac{1}{(1-r^2)^{3/2}}$, we get

$$ \dot{f}(r) = \gamma \underbrace{\int_{-\infty}^{\omega_i} \varphi(\alpha - \beta x_i) \varphi(x_i) dx_i}_{J} - \delta \underbrace{\int_{-\infty}^{\omega_i} x_i \varphi(\alpha - \beta x_i) \varphi(x_i) dx_i}_{I}. $$

Yet

$$ \begin{aligned} I &= - \int_{-\infty}^{\omega_i} \varphi(\alpha - \beta x_i) \dot{\varphi}(x_i) dx_i \\ &= - \left( \left[ \varphi(\alpha - \beta x_i) \varphi(x_i) \right]_{-\infty}^{\omega_i} + \beta \int_{-\infty}^{\omega_i} \dot{\varphi}(\alpha - \beta x_i) \varphi(x_i) dx_i \right) \\ &= -\varphi(\alpha - \beta \omega_i) \varphi(\omega_i) + \beta \int_{-\infty}^{\omega_i} (\alpha - \beta x_i) \varphi(\alpha - \beta x_i) \varphi(x_i) dx_i, \end{aligned} $$

hence

$$ I = -\varphi(\omega_i)\varphi(\alpha - \beta\omega_i) - \beta^2 I + \alpha\beta J $$

i.e.

$$ I = \frac{-\varphi(\omega_i)\varphi(\alpha - \beta\omega_i)}{1 + \beta^2} + \frac{\alpha\beta}{1 + \beta^2}J. $$

We thus have

$$ \begin{aligned} \dot{f}(r) &= \gamma J + \frac{\delta \varphi(\omega_i) \varphi(\alpha - \beta \omega_i)}{1 + \beta^2} - \frac{\alpha \beta}{1 + \beta^2} J \\ &= \frac{\delta \varphi(\omega_i) \varphi(\alpha - \beta \omega_i)}{1 + \beta^2} + \left( \gamma - \frac{\alpha \beta \delta}{1 + \beta^2} \right) J, \end{aligned} $$

but in the other hand $\frac{\delta}{1 + \beta^2} = \frac{1}{\sqrt{1 - r^2}}$ and

$$ \gamma - \frac{\alpha\beta\delta}{1 + \beta^2} = \frac{r\alpha}{1 - r^2} - \frac{\alpha r}{1 - r^2} = 0. $$

Thus $\dot{f}(r) = \frac{\varphi(\omega_i)\varphi(\alpha - \beta\omega_i)}{\sqrt{1 - r^2}} > 0$ and $f$ is strictly increasing. □
▷ Assumption **H5** is already assumed here.

▷ In order to verify assumptions **H6**, and **H7**, we denote once again for any $\theta \in \Theta$ and any $Z_0$

$$m_{\theta}(Z_0) = \log(f_{\theta}(Z_0)),$$

where

$$f_{\theta}(Z_0) = \int_{\mathbb{R}^k} \prod_{i=1}^{k} 1_{I_{Y_{i,0}} - \lambda_{i,0}(\gamma)}(x_i) \varphi_R(x) dx.$$

Let $\theta_\alpha$ and $\theta_\beta$ be two coordinates of the vector $\theta$. If the partial derivatives exist, we actually have

$$\frac{\partial m_{\theta}}{\partial \theta_{\alpha}}(Z_0) = \frac{\partial f_{\theta}}{\partial \theta_{\alpha}}(Z_0) \times \frac{1}{f_{\theta}(Z_0)}$$

and

$$\frac{\partial^2 m_{\theta}}{\partial \theta_{\beta} \partial \alpha}(Z_0) = \frac{\partial^2 f_{\theta}}{\partial \theta_{\beta} \partial \alpha}(Z_0) \times \frac{1}{f_{\theta}(Z_0)} - \frac{\partial f_{\theta}}{\partial \theta_{\beta}}(Z_0) \times \frac{\partial f_{\theta}}{\partial \theta_{\alpha}}(Z_0) \times \frac{1}{f_{\theta}(Z_0)^2}.$$

For future reasons, we will establish the stronger following result, which implies immediately the verification of assumption **H6** and **H7**.

**Lemma 7.5.** Let $q \in (1, +\infty)$, the application $\theta \mapsto m_{\theta}(Z_0)$ is of class $C^2$ on $\Theta$ and we have

$$\mathbb{E} \left( \sup_{\theta \in \Theta} \left| \frac{\partial m_{\theta}}{\partial \theta_{\alpha}}(Z_0) \right|^q \right) < +\infty \quad \text{and} \quad \mathbb{E} \left( \sup_{\theta \in \Theta} \left| \frac{\partial^2 m_{\theta}}{\partial \theta_{\beta} \partial \theta_{\alpha}}(Z_0) \right|^q \right) < +\infty.$$

The proof of Lemma 7.5 is not exactly straightforward, and requires the following preliminary result.

**Lemma 7.6.** The application $\theta \mapsto f_{\theta}(Z_0)$ is of class $C^2$ on $\Theta$.

Moreover, there exist two polynomials $P, Q$ with positive coefficients, such that for any $\theta \in \Theta$, we have

$$\left| \frac{\partial f_{\theta}}{\partial \theta_{\alpha}}(Z_0) \right| \le \int_{\mathbb{R}^k} \prod_{i=1}^{k} 1_{I_{Y_{i,0}} - \lambda_{i,0}(\gamma)}(x_i) [P(\|x\|_2)Q(\|X_{-1}\|_2)] \varphi_R(x) dx \quad (26)$$

and

$$\left| \frac{\partial^2 f_{\theta}}{\partial \theta_{\beta} \partial \theta_{\alpha}}(Z_0) \right| \le \int_{\mathbb{R}^k} \prod_{i=1}^{k} 1_{I_{Y_{i,0}} - \lambda_{i,0}(\gamma)}(x_i) [P(\|x\|_2)Q(\|X_{-1}\|_2)] \varphi_R(x) dx \quad (27)$$

*Proof of Lemma 7.6.* Due to the number of different cases to deal with, the entire proof of this lemma is quite repetitive and tedious. We only show here one specific case, but the other ones can be treated with the same arguments.

Assume $\theta_\alpha = B(i_0, j_0)$, we actually have

$$\begin{aligned} f_{\theta}(Z_0) &= \int_{\mathbb{R}^k} \prod_{i=1}^{k} 1_{I_{Y_{i,0}} - \lambda_{i,0}(\gamma)}(x_i) \varphi_R(x) dx \\ &= \int_{\mathbb{R}^{k-1}} \prod_{i \neq i_0} 1_{I_{Y_{i,0}} - \lambda_{i,0}(\gamma)}(x_i) \frac{\exp\left(-\frac{1}{2} \langle \overline{R^{-1}i_0} \bar{x}^{i_0}, \bar{x}^{i_0} \rangle\right)}{\sqrt{(2\pi)^k \det(R)}} \\ &\quad \times \int_{\mathbb{R}} 1_{I_{Y_{i_0,0}} - \lambda_{i_0,0}(\gamma)}(x_{i_0}) \exp\left(-\frac{1}{2} \overline{R^{-1}(i_0, i_0)x_{i_0}^2} - \sum_{i \neq i_0} \overline{R^{-1}(i_0, i)x_{i_0}x_i}\right) dx_{i_0} d\bar{x}^{i_0} \end{aligned}$$

where

* $\bar{x}^{i_0}$ denotes the vector $x$ deprived of the $i_0^{th}$ coordinate,
• $\overline{R^{-1}i_0}$ denotes the matrix $R^{-1}$ deprived of the $i_0^{th}$ line and $i_0^{th}$ column.

Denoting $g_{\theta,Z_0}^{(i_0)}(\bar{x}^{i_0})$ the integrand of this latter integral, we have if $X_{-1}(j_0) = 0$

$$ \frac{\partial}{\partial \theta_{\alpha}} \{ g_{\theta, Z_0}^{(i_0)}(\bar{x}^{i_0}) \} = 0 $$

and

$$ \frac{\partial}{\partial \theta_{\alpha}} \{ g_{\theta, Z_0}^{(i_0)}(\bar{x}^{i_0}) \} = (2Y_{i_0} - 1)X_{-1}(j_0) \prod_{i \neq i_0} 1_{I_{Y_{i_0,0}} - \lambda_{i_0,0}(\gamma)}(x_i) \varphi_R(x_{\lambda_{i_0,0}(\gamma)}) $$

otherwise, where $x_{\lambda_{i_0,0}(\gamma)} = (x_1, \dots, x_{i_0-1}, -\lambda_{i_0,0}(\gamma), x_{i_0+1}, \dots, x_k)$. Using Lemma 7.3, there exist positive constants $m_{\text{det}}, c_{20}$ such that for any $\theta \in \Theta$ and any $\bar{x}^{i_0} \in \mathbb{R}^{k-1}$

$$ \left| \frac{\partial}{\partial \theta_{\alpha}} \{ g_{\theta, Z_0}^{(i_0)}(\bar{x}^{i_0}) \} \right| \le \frac{\|X^{-1}\|_2}{\sqrt{(2\pi)^k m_{\text{det}}}} \exp \left( -\frac{c_2}{2} \|x_{\lambda_{i_0,0}(\gamma)}\|_2^2 \right) \\ \le \frac{\|X^{-1}\|_2}{\sqrt{(2\pi)^k m_{\text{det}}}} \exp \left( -\frac{c_2}{2} \|\bar{x}^{i_0}\|_2^2 \right) $$

which is integrable on $\mathbb{R}^{k-1}$. Therefore, the partial derivative $\frac{\partial f_{\theta}}{\partial \theta_{\alpha}}(Z_0)$ exists, and is either equal to zero or

$$ \frac{\partial f_{\theta}}{\partial \theta_{\alpha}}(Z_0) = (2Y_{i_0} - 1)X_{-1}(j_0) \int_{\mathbb{R}^{k-1}} \prod_{i \neq i_0} 1_{I_{Y_{i_0,0}-\lambda_{i_0,0}(\gamma)}}(x_i) \varphi_R(x_{\lambda_{i_0,0}(\gamma)}) d\bar{x}^{i_0}. $$

In order to obtain the inequality (26), we can use the fundamental theorem of calculus, and we obtain for any $\theta \in \Theta$

$$ \begin{aligned} \varphi_R(x_{\lambda_{i_0,0}(\gamma)}) &= \int_{-\infty}^{-\lambda_{i_0,0}(\gamma)} \frac{d}{dx_{i_0}} \{\varphi_R(x)\} dx_{i_0} \\ &= \left| \int_{I_{Y_{i_0,0}-\lambda_{i_0,0}(\gamma)}} \frac{d}{dx_{i_0}} \{\varphi_R(x)\} dx_{i_0} \right| \\ &\le \int_{I_{Y_{i_0,0}-\lambda_{i_0,0}(\gamma)}} |\langle R^{-1}(i_0, \cdot), x \rangle \varphi_R(x)| dx_{i_0} \\ &\le \int_{I_{Y_{i_0,0}-\lambda_{i_0,0}(\gamma)}} M_{\mathbb{R}^{-1}} \|x\|_2 \varphi_R(x) dx_{i_0} \end{aligned} $$

where $M_{\mathbb{R}^{-1}}$ is a positive constant, obtained by compactness of $\mathbb{R}$. Therefore, for any $\theta \in \Theta$

$$ \left| \frac{\partial f_{\theta}}{\partial \theta_{\alpha}}(Z_0) \right| \le M_{\mathbb{R}^{-1}} \int_{\mathbb{R}^k} \prod_{i=1}^k 1_{I_{Y_{i_0,0}-\lambda_{i_0,0}(\gamma)}}(x_i) \times \|x\|_2 \times \|X_{-1}\|_2 \times \varphi_R(x) dx, $$

which leads to the desired result.

Every other partial derivative, up to order 2, can be bounded the same way, and we obtain Lemma 7.6. $\square$

*Proof of Lemma 7.5.* Take any coordinates $\theta_{\alpha}, \theta_{\beta}$ of the vector $\theta$. By Lemma 7.6, the application $\theta \mapsto m_{\theta}(Z_0)$ is of class $C^2$ and we have

$$ \frac{\partial m_{\theta}}{\partial \theta_{\alpha}}(Z_0) = \frac{\partial f_{\theta}}{\partial \theta_{\alpha}}(Z_0) \times \frac{1}{f_{\theta}(Z_0)} $$

and

$$ \frac{\partial^2 m_{\theta}}{\partial \theta_{\beta} \partial_{\alpha}}(Z_0) = \frac{\partial^2 f_{\theta}}{\partial \theta_{\beta} \partial_{\alpha}}(Z_0) \times \frac{1}{f_{\theta}(Z_0)} - \frac{\partial f_{\theta}}{\partial \theta_{\beta}}(Z_0) \times \frac{\partial f_{\theta}}{\partial \theta_{\alpha}}(Z_0) \times \frac{1}{f_{\theta}(Z_0)^2}. $$
It is however possible to apply Hölder's inequality with respect to the measure of density $\varphi_R$ to inequality (26), and we obtain for any $\theta \in \Theta$

$$ \left| \frac{\partial f_{\theta}}{\partial \theta_{\alpha}}(Z_0) \right| \le C \int_{\mathbb{R}^k} \prod_{i=1}^k 1_{I_{Y_{i,0}} - \lambda_{i,0}(\gamma)}(x_i) [P(\|x\|_2)Q(\|X_{-1}\|_2)] \varphi_R(x) dx \\ \le C \left( \underbrace{\int_{\mathbb{R}^k} \prod_{i=1}^k 1_{I_{Y_{i,0}} - \lambda_{i,0}(\gamma)}(x_i) \varphi_R(x) dx}_{f_{\theta}(Z_0)} \right)^{\frac{1}{v}} \times \left( \underbrace{\int_{\mathbb{R}^k} [P(\|x\|_2)Q(\|X_{-1}\|_2)]^w \varphi_R(x) dx}_{I_{\theta}(X_{-1})} \right)^{\frac{1}{w}} $$

where $\frac{1}{v} + \frac{1}{w} = 1$, and $v, w > 1$ are to be determined later. Using Lemma 7.3, we have

$$ I_{\theta}(X_{-1}) \le Q(\|X_{-1}\|_2)^w \times \frac{1}{\sqrt{(2\pi)^k m_{\text{det}}}} \int_{\mathbb{R}^k} P(\|x\|_2)^w \exp\left(-\frac{c_2}{2}\|x\|_2^2\right) dx. $$

Since the Gaussian distribution admits moments of any order, one can find a positive constant $M_w$ depending only on $w$ such that for any $\theta \in \Theta$

$$ I_{\theta}(X_{-1})^{\frac{1}{w}} \le M_w Q(\|X_{-1}\|_2). $$

Thus, we have for any $\theta \in \Theta$

$$ \left| \frac{\partial m_{\theta}}{\partial \theta_{\alpha}}(Z_0) \right| \le CM_w Q(\|X_{-1}\|_2) \times \left( \int_{\mathbb{R}^k} \prod_{i=1}^k 1_{I_{Y_{i,0}} - \lambda_{i,0}(\gamma)}(x_i) \varphi_R(x) dx \right)^{\frac{1}{v}-1}. $$

Yet, by Lemma 7.3, there exist positive constants $M_{\text{det}}$ and $c_1$ such that for any $\theta \in \Theta$

$$ \begin{aligned} \int_{\mathbb{R}^k} \prod_{i=1}^k 1_{I_{Y_{i,0}} - \lambda_{i,0}(\gamma)}(x_i) \varphi_R(x) dx &\ge \frac{1}{\sqrt{M_{\text{det}}} \int_{\mathbb{R}^k} \prod_{i=1}^k 1_{I_{Y_{i,0}} - \lambda_{i,0}(\gamma)}(x_i) \exp\left(-\frac{c_1}{2}\|x\|_2^2\right) dx} \\ &= \frac{1}{\sqrt{c_1 M_{\text{det}}} \int_{\mathbb{R}^k} \prod_{i=1}^k 1_{I_{Y_{i,0}} - \sqrt{c_1} \lambda_{i,0}(\gamma)}(x_i) \exp\left(-\frac{\|x\|_2^2}{2}\right) dx} \\ &\ge \frac{1}{\sqrt{c_1 M_{\text{det}}} \prod_{i=1}^k \min\{\Phi(-\sqrt{c_1} \lambda_{i,0}(\gamma)), 1 - \Phi(-\sqrt{c_1} \lambda_{i,0}(\gamma))\}} \\ &\ge \frac{d}{\sqrt{c_1 M_{\text{det}}} \exp(-c_1 \|\lambda_0(\gamma)\|_2^2)} \end{aligned} $$

where $d$ is the positive constant from Lemma 7.8

$$ \ge \frac{d}{\sqrt{c_1 M_{\text{det}}} \exp(-c_1 (M_A + M_B \|X_{-1}\|_2)^2)} $$

for some positive constants $M_A$ and $M_B$, obtained by compactness of $\Theta$. We finally obtain that for any $\theta \in \Theta$

$$ \left| \frac{\partial m_{\theta}}{\partial \theta_{\alpha}}(Z_0) \right| \le CM_w Q(\|X_{-1}\|_2) \times \left( \frac{d}{\sqrt{c_1 M_{\text{det}}}} \right)^{\frac{1}{v}-1} \exp \left( c_1 \left( 1 - \frac{1}{v} \right) (M_A + M_B \|X_{-1}\|_2)^2 \right), $$

and for any $q > 1$

$$ \sup_{\theta \in \Theta} \left| \frac{\partial m_{\theta}}{\partial \theta_{\alpha}}(Z_0) \right|^q \le C^q M_w^q Q(\|X_{-1}\|_2)^q \times \left( \frac{d}{\sqrt{c_1 M_{\text{det}}}} \right)^{q(\frac{1}{v}-1)} \\ \times \exp \left( q c_1 \left( 1 - \frac{1}{v} \right) (M_A + M_B \|X_{-1}\|_2)^2 \right). $$
However, $v > 1$ is arbitrary, and we can choose it such that $qc_1(1 - \frac{1}{v}) M_B^2 < \kappa$, where $\kappa$ satisfies (17)
We thus have for any $q > 1$

$$
\mathbb{E} \left( \sup_{\theta \in \Theta} \left| \frac{\partial m_{\theta}}{\partial \theta_{\alpha}}(Z_0) \right|^q \right) < +\infty.
$$

Using inequality (27) from Lemma 7.6, we can show similarly that for any $q > 1$

$$
\mathbb{E} \left( \sup_{\theta \in \Theta} \left| \frac{\partial^2 m_{\theta}}{\partial \theta_{\alpha} \partial \theta_{\beta}}(Z_0) \right|^q \right) < +\infty.
$$

□

▷ Let us now check **H8**. According to Lemma 7.5, the application

$$
\theta \mapsto \mathbb{E}(m_{\theta}(Z_0))
$$

is twice differentiable. Thus, we can write the Taylor expansion for all $\theta \in \Theta$

$$
\mathbb{E}(m_{\theta}(Z_0)) = \mathbb{E}(m_{\theta_0}(Z_0)) + \mathbb{E}(\nabla m_{\theta_0}(Z_0)) \cdot (\theta - \theta_0) + (\theta - \theta_0)' \cdot \mathbb{E}(\nabla^2 m_{\theta_0}(Z_0)) \cdot (\theta - \theta_0) + o(\|\theta - \theta_0\|_2^2).
$$

Additionally, since $\theta_0$ is the unique maximum of the latter function, we actually have

$$
\mathbb{E}(m_{\theta}(Z_0)) = \mathbb{E}(m_{\theta_0}(Z_0)) + (\theta - \theta_0)' \cdot \mathbb{E}(\nabla^2 m_{\theta_0}(Z_0)) \cdot (\theta - \theta_0) + o(\|\theta - \theta_0\|_2^2),
$$

and the matrix $\mathbb{E}(\nabla^2 m_{\theta_0}(Z_0))$ is thus symmetric semi positive definite. It is then invertible if and only if it is positive definite. As usual for regular conditional likelihood estimators, one can show that non-invertibility of $M := \mathbb{E}(\nabla^2 m_{\theta_0}(Z_0))$ is equivalent to the existence of a linear relation between the partial derivatives of $\theta \mapsto f_{\theta}(Z_0)$ at point $\theta_0$. Let us define

$$
W_t = (1', Y'_{t-1}, \dots, Y'_{t-p}, X'_{t-1})',
$$

where $\mathbf{1}$ is the vector of $\mathbb{R}^k$ with all the coordinates equal to 1 and $D$ be the matrix of size $k \times (p+1)k+d$
such that $\lambda_t(\gamma) = DW_t$. Moreover if $D_1, \dots, D_k$ denote the row vectors of $D$, we have $\lambda_{i,t}(\gamma_0) = D'_i W_t$
for $1 \le i \le d$. When $M$ is non-invertible, one can find some vectors $\alpha_1, \dots, \alpha_k$ of $\mathbb{R}^{(p+1)k+d}$ and some
real numbers $\kappa_{m,l}$ for $1 \le m < l \le k$ such that a.s.

$$
0 = \sum_{i=1}^{k} (2Y_{i,1} - 1)\alpha'_{i}W_{0} \int_{\prod_{j \neq i} \{I_{Y_{j,1}} - D'_{j}W_{0}\}} \phi_{R} (x_{1}, \dots, x_{i-1}, B'_{i}W_{0}, x_{i+1}, \dots, x_{k}) \prod_{j \neq i} dx_{j} \\
+ \int_{\prod_{i=1}^{k} \{I_{Y_{i,1}} - Z'_{i}W_{0}\}} \sum_{1 \le m < l \le k} \kappa_{m,l} \frac{\partial \phi_{R}}{\partial R_{m,l}}(x)dx.
$$

Since $\mathbb{P}(Y_t = y | \mathcal{F}_{t-1})$ is positive a.s. for any value of $y \in \{0, 1\}^k$, we then get

$$
0 = \sum_{i=1}^{k} (2y_{i,1} - 1)\alpha'_{i}W_{0} \int_{\prod_{j \neq i} \{I_{y_{j,1}} - D'_{j}W_{0}\}} \phi_{R} (x_{1}, \dots, x_{i-1}, B'_{i}W_{0}, x_{i+1}, \dots, x_{k}) \prod_{j \neq i} dx_{j} \\
+ \int_{\prod_{i=1}^{k} \{I_{y_{i,1}} - Z'_{i}W_{0}\}} \sum_{1 \le m < l \le k} \kappa_{m,l} \frac{\partial \phi_{R}}{\partial R_{m,l}}(x)dx.
$$

almost surely. But summing these relations over $y_{j,1} \in \{0,1\}^k$ for $j \neq i$, we obtain

$$
\alpha'_{i} W_{0} = 0, \quad 1 \le i \le k, \text{ a.s.}
$$

As for the proof for consistency with the identification of $\theta_0$, one can show that $\alpha_i = 0$ for $1 \le i \le k$.
Summing now the previous relations over $y_{\ell,1}$ for $\ell \ne i, j$ for some $1 \le i < j \le k$, we get

$$
\kappa_{i,j} \int_{I_{y_{i,1}} - D'_{i} W_{0}} \int_{I_{y_{j,1}} - D'_{j} W_{0}} \frac{\partial \phi_{R}}{\partial R_{i,j}}(x_i, x_j) dx_i dx_j = 0 \text{ a.s.}
$$

However, the latter double integral cannot vanish because it is the derivative of a monotonic function
(it is equal to the function $f$ defined in Lemma 7.4 if we set $y_{i,1} = y_{j,1} = 0$, $r = R_{i,j}$ and $\omega_\ell = -D'_\ell W_0$
for $\ell = i, j$). Then $\kappa_{i,j} = 0$ and Assumption **H8** is satisfied.
▷ Finally, denoting once again

$$ \psi_{\theta}(y_t, Y_{t-1}^{-}, X_{t-1}) = \int_{\mathbb{R}^k} \prod_{i=1}^{k} \mathbb{1}_{I_{y_{i,t}}}(\lambda_{i,t}(\gamma) + e_i) \varphi_R(e) \mathrm{d}e $$

the density of the conditional distribution $\mathcal{L}_{\theta}(Y_t | \mathcal{F}_{t-1})$ with respect to the measure $\mu = \sum_{y \in \{0,1\}^k} \delta_y$, we have

$$ \nabla m_{\theta}(Z_t) = \frac{\nabla \psi_{\theta}(Y_t, Y_{t-1}^{-}, X_{t-1})}{\psi_{\theta}(Y_t, Y_{t-1}^{-}, X_{t-1})}. $$

Differentiating with respect to $\theta$ both sides of

$$ \int \psi_{\theta}(y_t, Y_{t-1}^{-}, X_{t-1}) \mu(dy_t) = 1, $$

we obtain for all $\theta \in \Theta$

$$ \int \nabla \psi_{\theta}(y_t, Y_{t-1}^{-}, X_{t-1}) \mu(dy_t) = 0. $$

Note that this differentiation is licit according to what we established in the previous points, since $\psi_{\theta}(Y_t, Y_{t-1}^{-}, X_{t-1}) = f_{\theta}(Z_t)$.

Thus, for any positive measurable function $h$ defined on $(\mathbb{R}^d)^N \times (\mathbb{R}^k)^N$, we have

$$ \begin{aligned} \mathbb{E} (h(X_{t-1}^{-}, \varepsilon_{t-1}^{-}) \nabla m_{\theta_0}(Z_t)) &= \iint \nabla m_{\theta_0}(y_t, \dots, y_{t-p}, x_{t-1}) \psi_{\theta_0}(y_t) \mu(dy_t) \\ & \quad h(x_{t-1}^{-}, e_{t-1}^{-}) \mathbb{P}_{(X_{t-1}^{-}, \varepsilon_{t-1}^{-})} (dx_{t-1}^{-}, de_{t-1}^{-}) \\ &= \iint \nabla \psi_{\theta_0}(y_t, y_{t-1}^{-}, x_{t-1}) \mu(dy_t) \\ & \quad h(x_{t-1}^{-}, e_{t-1}^{-}) \mathbb{P}_{(X_{t-1}^{-}, \varepsilon_{t-1}^{-})} (dx_{t-1}^{-}, de_{t-1}^{-}) \\ &= 0. \end{aligned} $$

Thus, $\mathbb{E}(\nabla m_{\theta_0}(Z_t) | \mathcal{F}_{t-1}) = 0$, i.e. Assumption **H9** is satisfied.

### 7.3.2 Proof for estimator $\hat{\gamma}$

The proof for estimator $\hat{\gamma}$ is quite similar to the one for $\hat{\theta}_{\mathcal{P}_L}$. Once again we have to check assumptions **H1** to **H9** from Theorem 7.9 and Theorem 7.10, for the same process $(Z_t)_{t \in \mathbb{Z}}$ and the family of applications $(m_{\gamma})_{\gamma \in \Gamma}$ defined by

$$ m_{\gamma}(Z_t) = \sum_{i=1}^{k} Y_{i,t} \log[\Phi(\lambda_{i,t}(\gamma))] + (1 - Y_{i,t}) \log[\Phi(-\lambda_{i,t}(\gamma))]. $$

▷ **H1** and **H2** are satisfied for the same reasons as before.

▷ The application $\gamma \mapsto m_{\gamma}(Z_0)$ is continuous almost surely by a simple composition of continuous functions.

The verification of the second point of **H3** relies once again on a compactness argument and the following inequality (see Lemma 7.8)

$$ \forall x \in \mathbb{R}, \quad 1 - \Phi(x) \geq de^{-x^2} \qquad (28) $$

for some constant $d > 0$.

Indeed, there exist positive constants $M_A$ and $M_B$ such that for all $\gamma \in \Gamma$ and all $i \in \{1, \dots, k\}$

$$ |\lambda_{i,0}(\gamma)| \leq M_A + M_B \|X_{-1}\|_2. $$
We thus have for all $\gamma \in \Gamma$

$$ egin{aligned} |m_{\gamma}(Z_0)| &= - \left\{ \sum_{i=1}^{k} Y_{i,0} \log[\Phi(\lambda_{i,0}(\gamma))] + (1 - Y_{i,0}) \log[\Phi(-\lambda_{i,0}(\gamma))] \right\} \\ &\le -2k \log[\Phi(-M_A - M_B \|X_{-1}\|_2)] \\ &= -2k \log[1 - \Phi(M_A + M_B \|X_{-1}\|_2)]. \end{aligned} $$

Then, from inequality (28), we obtain

$$ egin{aligned} |m_{\gamma}(Z_0)| &\le -2k \log(d \exp(-(M_A + M_B \|X_{-1}\|_2)^2)) \\ &\le -2k \log(d) + 2k (M_A + M_B \|X_{-1}\|_2)^2. \end{aligned} $$

Hence

$$ \mathbb{E} \left( \sup_{\gamma \in \Gamma} |m_{\gamma}(Z_0)| \right) < +\infty $$

and assumption **H3** is satisfied.

▷ Using the same notations as before, we have

$$ \mathbb{E}(m_{\gamma}(Z_0)) - \mathbb{E}(m_{\gamma_0}(Z_0)) = \sum_{i=1}^{k} \mathbb{E} \left( \log \left( \frac{\psi_{\gamma}(Y_{i,0}, Y_{-1}^{-}, X_{-1})}{\psi_{\gamma_0}(Y_{i,0}, Y_{-1}^{-}, X_{-1})} \right) \right), $$

where by assumption **B2**, $\psi_{\gamma}(y_{i,0}, Y_{-1}, X_{-1}) = \Phi(\lambda_{i,0}(\gamma))^{y_{i,0}} \Phi(-\lambda_{i,0}(\gamma))^{1-y_{i,0}}$ is the density of the conditional distribution $\mathcal{L}_{\gamma}(Y_{i,0} | Y_{-1}^{-}, X_{-1})$ with respect to the measure $\mu = \delta_0 + \delta_1$ on $\{0, 1\}$.

The same reasoning as before ensures that we have $\mathbb{E}(m_{\gamma}(Z_0)) \le \mathbb{E}(m_{\gamma_0}(Z_0))$ with equality if and only if $\gamma = \gamma_0$, so assumption **H4** is verified.

▷ Assumption **H5** is already assumed here.

▷ A simple computation gives for all $(i_1, i_2) \in \{1, \dots, k\}^2$ and all $l \in \{1, \dots, p\}$

$$ \frac{\partial m_{\gamma}(Z_0)}{\partial A_l(i_1, i_2)} = Y_{-l, i_2} \left( Y_{i_1, 0} \frac{\varphi(\lambda_{i_1, 0}(\gamma))}{\Phi(\lambda_{i_1, 0}(\gamma))} + (Y_{i_1, 0} - 1) \frac{\varphi(-\lambda_{i_1, 0}(\gamma))}{\Phi(-\lambda_{i_1, 0}(\gamma))} \right), $$

which leads to

$$ \left| \frac{\partial m_{\gamma}(Z_0)}{\partial A_l(i_1, i_2)} \right|^2 = Y_{-l, i_2} \times \frac{\varphi(\lambda_{i_1, 0}(\gamma))^2}{\Phi(\lambda_{i_1, 0}(\gamma))^2} \mathbb{1}_{\lambda_{i_1, 0}(\gamma)>0} + Y_{-l, i_2} \times \frac{\varphi(-\lambda_{i_1, 0}(\gamma))^2}{\Phi(-\lambda_{i_1, 0}(\gamma))^2} \mathbb{1}_{\lambda_{i_1, 0}(\gamma)\le 0}. $$

Yet, for any $x \ge 0$, we have $\Phi(x) \ge \frac{1}{2}$, thus $\frac{\varphi(x)^2}{\Phi(x)^2} \le 4$, and it follows

$$ \sup_{\gamma \in \Gamma} \left| \frac{\partial m_{\gamma}(Z_0)}{\partial A_l(i_1, i_2)} \right|^2 \le 4. $$

Similarly, we have for any $i_1 \in \{1, \dots, k\}$ and $i_2 \in \{1, \dots, d\}$

$$ \frac{\partial m_{\gamma}(Z_0)}{\partial B(i_1, i_2)} = X_{i_2, -1} \left( Y_{i_1, 0} \frac{\varphi(\lambda_{i_1, 0}(\gamma))}{\Phi(\lambda_{i_1, 0}(\gamma))} + (Y_{i_1, 0} - 1) \frac{\varphi(-\lambda_{i_1, 0}(\gamma))}{\Phi(-\lambda_{i_1, 0}(\gamma))} \right), \quad (29) $$

and we obtain

$$ \sup_{\gamma \in \Gamma} \left| \frac{\partial m_{\gamma}(Z_0)}{\partial B(i_1, i_2)} \right|^2 \le 4 \|X_{-1}\|_2^2. $$

Since $\mathbb{E}(\|X_{-1}\|_2^2) < +\infty$, we have

$$ \mathbb{E} \left( \sup_{\gamma \in \Gamma} \| \nabla m_{\gamma}(Z_0) \|_2^2 \right) < +\infty, $$

and **H6** is satisfied.
▷ Setting for any $x \in \mathbb{R}$, $g(x) = \frac{\varphi(x)}{\Phi(x)}$, and using equality (29), we have for any $i_1, i_2, i_3, i_4$ in $\{1, \dots, k\}$
and any $l, l'$ in $\{1, \dots, p\}$

$$
\begin{aligned}
\frac{\partial^2 m_\gamma(Z_0)}{\partial A_{l'}(i_3, i_4) \partial A_l(i_1, i_2)} &= Y_{-l, i_2} \left( Y_{i_1, 0} \frac{\partial g(\lambda_{i_1, 0}(\gamma))}{\partial A_{l'}(i_3, i_4)} + (Y_{i_1, 0} - 1) \frac{\partial g(-\lambda_{i_1, 0}(\gamma))}{\partial A_{l'}(i_3, i_4)} \right) \\
&= \begin{cases} Y_{-l, i_2} Y_{-l', i_4} (Y_{i_1, 0} \cdot g'(\lambda_{i_1, 0}(\gamma)) - (Y_{i_1, 0} - 1) \cdot g'(\lambda_{i_1, 0}(\gamma))) & \text{if } i_3 = i_1 \\ 0 & \text{otherwise.} \end{cases}
\end{aligned}
$$

Yet for any $x \ge 0$

$$
|g'(x)| = \frac{e^{-\frac{x^2}{2}}}{2\pi} \times \frac{x\sqrt{2\pi}\Phi(x) + e^{-\frac{x^2}{2}}}{\Phi(x)^2} \le \frac{2}{\pi} \times (1 + x\sqrt{2\pi}),
$$

thus we obtain

$$
\left| \frac{\partial^2 m_\gamma(Z_0)}{\partial A_{l'}(i_3, i_4) \partial A_l(i_1, i_2)} \right| \le \frac{2}{\pi} \left( 1 + |\lambda_{i_1, 0}(\gamma)| \sqrt{2\pi} \right).
$$

By a compactness argument, there exist two positive constants $M_A$ and $M_B$ such that for any $\gamma \in \Gamma$

$$
\left| \frac{\partial^2 m_\gamma(Z_0)}{\partial A_{l'}(i_3, i_4) \partial A_l(i_1, i_2)} \right| \le \frac{2}{\pi} \left( 1 + \sqrt{2\pi} (M_A + M_B \|X_{-1}\|_2) \right),
$$

and since $\mathbb{E}(\|X_{-1}\|_2) < +\infty$, we obtain

$$
\mathbb{E} \left( \sup_{\gamma \in \Gamma} \left| \frac{\partial^2 m_{\gamma}(Z_0)}{\partial A_{l'}(i_3, i_4) \partial A_l(i_1, i_2)} \right| \right) < +\infty.
$$

Similarly, we have for any $i_1, i_2, i_3$ in $\{1, \dots, k\}$, any $i_4 \in \{1, \dots, d\}$, and any $l \in \{1, \dots, p\}$

$$
\frac{\partial^2 m_\gamma(Z_0)}{\partial B(i_3, i_4) \partial A(i_1, i_2)} = \begin{cases} Y_{-l, i_2} X_{i_4, -1} (Y_{i_1, 0} \cdot g'(\lambda_{i_1, 0}(\gamma)) - (Y_{i_1, 0} - 1) \cdot g'(-\lambda_{i_1, 0}(\gamma))) & \text{if } i_1 = i_3 \\ 0 & \text{otherwise.} \end{cases}
$$

With the same compactness argument as before we obtain that for any $\gamma \in \Gamma$

$$
\begin{aligned}
\left| \frac{\partial^2 m_\gamma(Z_0)}{\partial B(i_3, i_4) \partial A(i_1, i_2)} \right| &\le \|X_{-1}\|_2 \cdot |g'(\lambda_{i_1, 0}(\gamma))| \\
&\le \frac{2\|X_{-1}\|_2}{\pi} (1 + |\lambda_{i_1, 0}(\gamma)| \cdot \sqrt{2\pi}) \\
&\le \frac{2\|X_{-1}\|_2}{\pi} (1 + \sqrt{2\pi}(M_A + M_B \|X_{-1}\|_2))
\end{aligned}
$$

for some positive constants $M_A$ and $M_B$.

Since $\mathbb{E}(\|X_{-1}\|_2^2) < +\infty$, we then obtain

$$
\mathbb{E} \left( \sup_{\gamma \in \Gamma} \left| \frac{\partial^2 m_{\gamma}(Z_0)}{\partial B(i_3, i_4) \partial A(i_1, i_2)} \right| \right) < +\infty.
$$

Finally, we have for any $i_1, i_3$ in $\{1, \dots, k\}$ and any $i_2, i_4$ in $\{1, \dots, d\}$

$$
\frac{\partial^2 m_\gamma(Z_0)}{\partial B(i_3, i_4) \partial B(i_1, i_2)} = \begin{cases} X_{i_4, -1} X_{i_2, -1} (Y_{i_1, 0} \cdot g'(\lambda_{i_1, 0}(\gamma)) - (Y_{i_1, 0} - 1) \cdot g'(-\lambda_{i_1, 0}(\gamma))) & \text{if } i_3 = i_1 \\ 0 & \text{otherwise.} \end{cases}
$$

Yet, for any $x \ge 0$, we have $\Phi(x) \ge \frac{e^{-\frac{x^2}{2}}}{2}$, thus

$$
\forall x \ge 0, |g'(x)| = \left| -\frac{e^{-\frac{x^2}{2}}}{2\pi} \times \frac{x\sqrt{2\pi}\Phi(x) + e^{-\frac{x^2}{2}}}{\Phi(x)^2} \right| \le \frac{2}{\pi} + \sqrt{\frac{2}{\pi}}xe^{-\frac{x^2}{2}}.
$$
Furthermore, since $xe^{-\frac{x^2}{2}} \xrightarrow{x \to +\infty} 0$, there exists $x_0 > 0$ such that for all $x \ge x_0$, $xe^{-\frac{x^2}{2}} \le 1$. Hence, for all $\gamma \in \Gamma$ we have

$$ \left| \frac{\partial^2 m_\gamma(Z_0)}{\partial B(i_3, i_4) \partial B(i_1, i_2)} \right| \le \|X_{-1}\|_2^2 \cdot \left( \frac{2}{\pi} + \sqrt{\frac{2}{\pi}}(x_0 + 1) \right). $$

Since $\mathbb{E}(\|X_{-1}\|_2^2) < +\infty$, we then have

$$ \mathbb{E} \left( \sup_{\gamma \in \Gamma} \left| \frac{\partial^2 m_\gamma(Z_0)}{\partial B(i_3, i_4) \partial B(i_1, i_2)} \right| \right) < +\infty, $$

and it yields

$$ \mathbb{E} \left( \sup_{\gamma \in \Gamma} \| \nabla^2 m_\gamma(Z_0) \|_2 \right) < +\infty, $$

so assumption **H7** is satisfied.

▷ Assumption **H8** can be verified the exact same way as in the proof for estimator $\hat{\theta}_{\mathcal{P}\mathcal{L}}$.

▷ The verification of **H9** follows the same principle as in subsection 6.3.1, yet we have to verify that the application $\int \psi_\gamma(y_{i,t}, Y_{t-1}, X_{t-1})\mu(dy_{i,t})$ is differentiable with respect to $\gamma$, where

$$ \psi_\gamma(y_{i,t}, Y_{t-1}^{-}, X_{t-1}) = \Phi(\lambda_{i,t}(\gamma))^{y_{i,t}} \Phi(-\lambda_{i,t}(\gamma))^{1-y_{i,t}} $$

is the density of the conditional distribution $\mathcal{L}_\gamma(Y_{i,t} | Y_{t-1}^{-}, X_{t-1})$ with respect to the measure $\mu = \delta_0 + \delta_1$ on $\{0, 1\}$.

To do so, note that for any $l \in \{1, \dots, p\}$, and any $i_1, i_2$ in $\{1, \dots, k\}$, we have

$$ \frac{\partial \psi_{\gamma}(y_{i,t}, Y_{t-1}^{-}, X_{t-1})}{\partial A_l(i_1, i_2)} = (Y_{i_2, t-l} \varphi(\lambda_{i,t}(\gamma))) \mathbb{1}_{\lambda_{i,t}(\gamma)>0} - (Y_{i_2, t-l} \varphi(-\lambda_{i,t}(\gamma))) \mathbb{1}_{\lambda_{i,t}(\gamma)\le 0} $$

thus for any $\gamma \in \Gamma$

$$ \left| \frac{\partial \psi_{\gamma}(y_{i,t}, Y_{t-1}^{-}, X_{t-1})}{\partial A_l(i_1, i_2)} \right| \le 2. $$

Similarly, we have for any $i_1 \in \{1, \dots, k\}$ and $i_2 \in \{1, \dots, d\}$

$$ \frac{\partial \psi_{\gamma}(y_{i,t}, Y_{t-1}^{-}, X_{t-1})}{\partial B(i_1, i_2)} = (X_{i_2, t-1} \varphi(\lambda_{i,t}(\gamma))) \mathbb{1}_{\lambda_{i,t}(\gamma)>0} - (X_{i_2, t-1} \varphi(-\lambda_{i,t}(\gamma))) \mathbb{1}_{\lambda_{i,t}(\gamma)\le 0}, $$

thus for any $\gamma \in \Gamma$

$$ \left| \frac{\partial \psi_{\gamma}(y_{i,t}, Y_{t-1}^{-}, X_{t-1})}{\partial B(i_1, i_2)} \right| \le 2 \|X_{t-1}\|_2. $$

The differentiation of $\gamma \mapsto \int \psi_\gamma(y_{i,t}, Y_{t-1}^{-}, X_{t-1})\mu(dy_{i,t})$ is thus licit.

### 7.3.3 Proof for estimator $\hat{r}$

For any couple $(i, j) \in \{1, \dots, k\}$ with $i < j$, $r_{i,j}$ belongs to a compact set $S_{i,j}$ of $[-1, 1]$, and setting for $s \in S_{i,j}$

$$ \ell_{i,j}(\gamma, s, Z_t) = \log \left( \int_{I_{Y_{i,t}-\lambda_{i,t}(\gamma)}} \Phi \left[ (2Y_{i,t} - 1) \frac{\lambda_{j,t}(\gamma) + sx_i}{\sqrt{1-s^2}} \right] \varphi(x_i) dx_i \right), $$

we actually have

$$ \hat{r}_{i,j} = \underset{s \in S_{i,j}}{\operatorname{argmax}} \sum_{t=p+1}^{T} \ell_{i,j}(\hat{\gamma}_t, s, Z_t). $$

We first deal with the consistency of the estimator $\hat{r}_{i,j}$, $i < j$. To this end, we check that assumptions **I1** to **I3** of Proposition 7.11
▷ Assumption **I1** is satisfied since $(Z_t)_{t \in \mathbb{Z}}$ is ergodic.

▷ Using the continuity of $\gamma \mapsto \lambda_{i,t}(\gamma)$, $\Phi$, the boundedness of $\Phi$ and Lebesgue's dominated convergence theorem, one can deduce that the mapping

$$ (\gamma, s) \mapsto \int_{I_{Y_{i,t}-\lambda_{i,t}(\gamma)}} \Phi \left[ (2Y_{i,t} - 1) \frac{\lambda_{j,t}(\gamma) - sx_i}{\sqrt{1-s^2}} \right] \varphi(x_i) dx_i $$

is almost surely continuous.

▷ We then show that

$$ \mathbb{E} \left( \sup_{\gamma \in \Gamma, s \in S_{i,j}} |\ell_{i,j}(\gamma, s, Z_0)| \right) < +\infty. \quad (30) $$

By Lemma 7.8, there exists a positive number $d$ such that

$$ \forall x \in \mathbb{R}, \Phi(x) \ge d \exp(-x^2). $$

Thus, for any $\gamma \in \Gamma$ and any $s \in S_{i,j}$

$$ \Phi \left[ (2Y_{i,0} - 1) \frac{\lambda_{j,0}(\gamma) + sx_i}{\sqrt{1-s^2}} \right] \ge d \exp \left( - \frac{(\lambda_{j,0}(\gamma) + sx_i)^2}{1-s^2} \right) \ge d \exp \left( - \frac{2\lambda_{j,0}(\gamma)^2 + 2s^2x_i^2}{1-s^2} \right). $$

Setting $h(s) = \frac{2s^2}{1-s^2} + \frac{1}{2}$, we then get

$$ \begin{aligned} \ell_{i,j}(\gamma, s, Z_0) &\ge \log(d) - \frac{2\lambda_{j,0}(\gamma)^2}{1-s^2} + \log \left( \int_{I_{Y_{i,0}-\lambda_{i,0}(\gamma)}} \exp \left( -\frac{2s^2x_i^2}{1-s^2} \right) \varphi(x_i) dx_i \right) \\ &= \log(d) - \frac{2\lambda_{j,0}(\gamma)^2}{1-s^2} - \frac{1}{2}\log(2h(s)) + \log \left( \int_{I_{Y_{i,0}-\sqrt{2h(s)}\lambda_{i,0}(\gamma)}} \varphi(x_i) dx_i \right) \\ &\ge \log(d) - \frac{2\lambda_{j,0}(\gamma)^2}{1-s^2} - \frac{1}{2}\log(2h(s)) \\ &\phantom{\ge} + \log \left( \min \left\{ \Phi(-\sqrt{2h(s)}\lambda_{i,0}(\gamma)), 1 - \Phi(-\sqrt{2h(s)}\lambda_{i,0}(\gamma)) \right\} \right) \\ &\ge 2\log(d) - \frac{2\lambda_{j,0}(\gamma)^2}{1-s^2} - \frac{1}{2}\log(2h(s)) - 2h(s)\lambda_{i,0}(\gamma)^2. \end{aligned} $$

By a compactness argument, there exist $M_1, M_2 > 0$ such that

$$ \begin{aligned} \sup_{\gamma \in \Gamma, s \in S_{i,j}} |\ell_{i,j}(\gamma, s, Z_0)| &= \sup_{\gamma \in \Gamma, s \in S_{i,j}} \{-\ell_{i,j}(\gamma, s, Z_0)\} \\ &\le M_1 + M_2 \|X_{-1}\|_2^2 \end{aligned} $$

Since $\mathbb{E}(\|X_{-1}\|_2^2) < +\infty$ using our assumptions, we get (30), and assumption **I2** is satisfied.

▷ One can finally show the exact same way as it was done in subsections 6.3.1 and 6.3.2, that for any $r_{i,j} \in S_{i,j}$

$$ \mathbb{E}(\ell_{i,j}(\gamma_0, r_{i,j}, Z_0)) \le \mathbb{E}(\ell_{i,j}(\gamma, r_{0,i,j}, Z_0)). $$

Furthermore, the equality case leads to the equality

$$ \int_{-\infty}^{-\lambda_{i,0}(\gamma_0)} \Phi \left( \frac{-\lambda_{j,0}(\gamma_0) - r_{i,j}x}{\sqrt{1-r_{i,j}^2}} \right) \varphi(x) dx = \int_{-\infty}^{-\lambda_{i,0}(\gamma_0)} \Phi \left( \frac{-\lambda_{j,0}(\gamma_0) - r_{0,i,j}x}{\sqrt{1-r_{0,i,j}^2}} \right) \varphi(x) dx $$

which is the same as (25). We thus have $r_{i,j} = r_{0,i,j}$ in this case, and **I3** is satisfied.
From Proposition 7.11, we then deduce that

$$ \hat{r} \xrightarrow[T \to +\infty]{} r_0 \quad \text{a.s.} $$

The stochastic expansion of $\sqrt{T-p}(\hat{r}-r_0)$ can be established using Proposition 7.12. Note that

$$ \hat{r} = \underset{r \in \mathcal{R}}{\operatorname{argmax}} \sum_{p+1 \le t \le T} \ell_t(\hat{\gamma}, r). $$

where $\ell_t$ is defined in (10).

In what follows, we check that conditions **I4** to **I8** of this latter proposition are satisfied.

▷ Condition **I4** is satisfied by assumption.

▷ We then check assumption **I6** and prove that

$$ \mathbb{E} \left( \sup_{\theta \in \Theta} \| \nabla_{1,2} \ell_0(\theta) \| \right) < +\infty \quad \text{and} \quad \mathbb{E} \left( \sup_{\theta \in \Theta} \| \nabla_2^2 \ell_0(\theta) \| \right) < +\infty. \quad (33) $$

Using the fact that $\ell_0(\theta)$ is a sum of functions involving different parameters, we fix $i$ and $j$ with $1 \le i < j \le k$ and set for any $\gamma \in \Gamma$ and $s \in S_{i,j}$

$$ \tilde{\ell}_{i,j}(\gamma, s, Z_0) = \int_{I_{Y_{i,0}} - \lambda_{i,0}(\gamma)} \Phi(g_{j,x}(\gamma, s, Z_0)) \varphi(x) dx, $$

where $g_{j,x}(\gamma, s, Z_0) = (2Y_{j,0} - 1) \frac{\lambda_{j,0}(\gamma) + sx}{\sqrt{1-s^2}}$. This way, we actually have

$$ \ell_0(\theta) = \sum_{1 \le i < j \le k} \log (\tilde{\ell}_{i,j}(\gamma, r_{i,j}, Z_0)). $$

Next we denote by $\gamma_\alpha$ one of the coordinates of the vector $\gamma$. Note that $\gamma_\alpha$ appears in the expression of either $\lambda_{i,0}(\gamma)$ or $\lambda_{j,0}(\gamma)$ but not in both. First, computing the different partial derivatives, one can find a positive constant $C$ such that

$$ \left| \frac{\partial g_{j,x}}{\partial s}(\gamma, s, Z_0) \right| \le C (1 + |x| + \|X_{-1}\|_2), \quad \left| \frac{\partial g_{j,x}}{\partial \gamma_\alpha}(\gamma, s, Z_0) \right| \le C (1 + \|X_{-1}\|_2), $$

and

$$ \left| \frac{\partial^2 g_{j,x}}{\partial s \partial \gamma_\alpha}(\gamma, s, Z_0) \right| \le C (1 + \|X_{-1}\|_2), \quad \left| \frac{\partial^2 g_{j,x}}{\partial s^2}(\gamma, s, Z_0) \right| \le C (1 + |x| + \|X_{-1}\|_2). $$

We next compute the first order partial derivatives of $\tilde{\ell}_{i,j}(\gamma, s, Z_0)$. More precisely, one can apply the theorem for the differentiability of a function defined by an integral to get

$$ \frac{\partial \tilde{\ell}_{i,j}}{\partial s}(\gamma, s, Z_0) = \int_{I_{Y_{i,0}} - \lambda_{i,0}(\gamma)} \varphi(g_{j,x}(\gamma, s, Z_0)) \varphi(x) \frac{\partial g_{j,x}}{\partial s}(\gamma, s, Z_0) dx. $$

Moreover, if $\gamma_\alpha$ is a coordinate that appears in the expression of $\lambda_{i,0}(\gamma)$

$$ \frac{\partial \tilde{\ell}_{i,j}}{\partial \gamma_\alpha}(\gamma, s, Z_0) = (2Y_{i,0} - 1) \frac{\partial \lambda_{i,0}}{\partial \gamma_\alpha}(\gamma) \times \Phi(g_{j,x}(\gamma, s, Z_0)) \varphi(-\lambda_{i,0}(\gamma)). $$

If $\gamma_\alpha$ is a coordinate that appears in the expression of $\lambda_{j,0}(\gamma)$, we have

$$ \frac{\partial \tilde{\ell}_{i,j}}{\partial \gamma_\alpha}(\gamma, s, Z_0) = \int_{I_{Y_{i,0}} - \lambda_{i,0}(\gamma)} \varphi(g_{j,x}(\gamma, s, Z_0)) \varphi(x) \frac{\partial g_{j,x}}{\partial \gamma_\alpha}(\gamma, s, Z_0) dx. $$
For the second order partial derivatives of $\tilde{\ell}_{i,j}(\gamma, s, Z_0)$, we have

$$ \frac{\partial^2 \tilde{\ell}_{i,j}(\gamma, s, Z_0)}{\partial s^2} = - \int_{I_{Y_{i,0}-\lambda_{i,0}(\gamma)}} g_{j,x}(\gamma, s, Z_0) \varphi(g_{j,x}(\gamma, s, Z_0)) \varphi(x) \left| \frac{\partial g_{j,x}}{\partial s}(\gamma, s, Z_0) \right|^2 dx \\ + \int_{I_{Y_{i,0}-\lambda_{i,0}(\gamma)}} \varphi(g_{j,x}(\gamma, s, Z_0)) \varphi(x) \frac{\partial^2 g_{j,x}}{\partial s^2}(\gamma, s, Z_0) dx. $$

Moreover, if $\gamma_\alpha$ is a coordinate of $\gamma$ that appears in the expression of $\lambda_{i,0}(\gamma)$,

$$ \frac{\partial^2 \tilde{\ell}_{i,j}(\gamma, s, Z_0)}{\partial s \partial \gamma_\alpha} = (2Y_{i,0} - 1) \varphi(g_{j,-\lambda_{i,0}(\gamma)}(\gamma, s, Z_0)) \varphi(-\lambda_{i,0}(\gamma)) \frac{\partial g_{j,-\lambda_{i,0}(\gamma)}}{\partial s}(\gamma, s, Z_0) \frac{\partial \lambda_{i,0}}{\partial \gamma_\alpha}(\gamma). $$

Note that we can always modify this derivative using the following integral expression that can be obtained from the fundamental theorem of calculus

$$ (2Y_{i,0} - 1) \varphi(g_{j,-\lambda_{i,0}(\gamma)}(\gamma, s, Z_0)) \varphi(-\lambda_{i,0}(\gamma)) = \int_{I_{Y_{i,0}-\lambda_{i,0}(\gamma)}} \varphi(g_{j,x}(\gamma, s, Z_0)) \varphi(x) \times \left( \frac{(2Y_{j,0}-1)s}{\sqrt{1-s^2}} g_{j,x}(\gamma, s, Z_0) + x \right) dx. $$

Thus, if $\gamma_\alpha$ is a coordinate of $\gamma$ that appears in the expression of $\lambda_{i,0}(\gamma)$

$$ \frac{\partial^2 \tilde{\ell}_{i,j}(\gamma, s, Z_0)}{\partial s \partial \gamma_\alpha} = \int_{I_{Y_{i,0}-\lambda_{i,0}(\gamma)}} \varphi(g_{j,x}(\gamma, s, Z_0)) \varphi(x) \left( \frac{(2Y_{j,0}-1)s}{\sqrt{1-s^2}} g_{j,x}(\gamma, s, Z_0) + x \right) \times \frac{\partial g_{j,-\lambda_{i,0}(\gamma)}}{\partial s}(\gamma, s, Z_0) \frac{\partial \lambda_{i,0}}{\partial \gamma_\alpha}(\gamma) dx $$

On the other hand, when $\gamma_\alpha$ is a coordinate of $\gamma$ that appears in the expression of $\lambda_{j,0}(\gamma)$, we have the expression

$$ \frac{\partial^2 \tilde{\ell}_{i,j}(\gamma, s, Z_0)}{\partial s \partial \gamma_\alpha} = \int_{I_{Y_{i,0}-\lambda_{i,0}(\gamma)}} \varphi(g_{j,x}(\gamma, s, Z_0)) \varphi(x) \frac{\partial^2 g_{j,x}}{\partial s \partial \gamma_\alpha}(\gamma, s, Z_0) dx \\ - \int_{I_{Y_{i,0}-\lambda_{i,0}(\gamma)}} \varphi(g_{j,x}(\gamma, s, Z_0)) \varphi(x) g_{j,x}(\gamma, s, Z_0) \frac{\partial g_{j,x}}{\partial s}(\gamma, s, Z_0) \\ \times \frac{\partial g_{j,x}}{\partial \gamma_\alpha}(\gamma, s, Z_0) dx. $$

Consequently, we can write the following partial derivatives

$$ \frac{\partial^2}{\partial s^2} \log (\tilde{\ell}_{i,j}(\gamma, s, Z_0)) = \underbrace{\frac{\partial^2 \tilde{\ell}_{i,j}}{\partial s^2}(\gamma, s, Z_0) \times \frac{1}{\tilde{\ell}_{i,j}(\gamma, s, Z_0)}}_{I_{i,j}^{(1)}(\gamma,s,Z_0)} - \underbrace{\left(\frac{\partial \tilde{\ell}_{i,j}}{\partial s}(\gamma, s, Z_0) \times \frac{1}{\tilde{\ell}_{i,j}(\gamma, s, Z_0)}\right)^2}_{I_{i,j}^{(2)}(\gamma,s,Z_0)} $$

and

$$ \frac{\partial^2}{\partial s \partial \gamma_\alpha} \log (\tilde{\ell}_{i,j}(\gamma, s, Z_0)) = \underbrace{\frac{\partial^2 \tilde{\ell}_{i,j}}{\partial s \partial \gamma_\alpha}(\gamma, s, Z_0) \times \frac{1}{\tilde{\ell}_{i,j}(\gamma, s, Z_0)}}_{I_{i,j}^{(3)}(\gamma,s,Z_0)} - \underbrace{\frac{\partial \tilde{\ell}_{i,j}}{\partial \gamma_\alpha}(\gamma, s, Z_0) \times \frac{\partial \tilde{\ell}_{i,j}}{\partial s}(\gamma, s, Z_0) \times \frac{1}{\tilde{\ell}_{i,j}(\gamma, s, Z_0)^2}}_{I_{i,j}^{(4)}(\gamma,s,Z_0)}. $$
Our aim is now to bound uniformly on $\Theta$ each quantity $I_{i,j}^{(k)}(\gamma, s, Z_0)$, $1 \le k \le 4$, by a random variable that is integrable.

Using Lemma 7.8 and the upper-bounds for $g_{j,x}$ and its partial derivatives, we can show that the partial derivative (34) can be bounded, up to a constant, by an integral $I$ defined by

$$I = \int_{I_{Y_{i,0}} - \lambda_{i,0}(\gamma)} \Phi(g_{j,x}(\gamma, s, Z_0)) \varphi(x) (1 + |x|^q + \|X_{-1}\|_2^q) dx$$

for some positive integer $q$. Now we apply the Hölder inequality, with $v^{-1} + w^{-1} = 1$ and $v > 1$ being determined latter. We have

$$I \le (\tilde{\ell}_{i,j}(\gamma, s, Z_0))^{1/v} \left( \int_{\mathbb{R}} \varphi(x) (1 + |x|^q)^w dx \right)^{1/w} + \|X_{-1}\|_2^q \tilde{\ell}_{i,j}(\gamma, s, Z_0).$$

Using Lemma 7.7, there also exist $C_1, C_2 > 0$ such that

$$\inf_{\theta \in \Theta} \tilde{\ell}_{i,j}(\gamma, s, Z_0) \ge C_1 \exp(-C_2 \|X_{-1}\|_2^2).$$

By choosing $v$ sufficiently close to one and using our exponential moment assumption, we deduce that $I_{i,j}^{(1)}(\gamma, s, Z_0)$ is bounded uniformly on $\Theta$ by an integrable random variable.

One can show the same way that each quantity $I_{i,j}^{(k)}(\gamma, s, Z_0)$, $2 \le k \le 3$, is also bounded uniformly on $\Theta$ by an integrable random variable.

We then easily deduce (33), and assumption **I6** is satisfied.

▷ The condition

$$\mathbb{E} (\|\nabla_2 \ell_0(\theta_0)\|^2) < +\infty$$

can be obtained using similar arguments as for getting (33), so **I5** is satisfied.

▷ We then prove that $L_{2,2} = \mathbb{E} \left( \frac{\partial^2 \ell_0(\theta_0)}{\partial r^2} \right)$ is invertible.

Since the different terms in the sum $\ell_0(\theta_0)$ depend on different parameters, $L_{2,2}$ contains only its diagonal terms

$$\mathbb{E} \left( \frac{\partial^2 \ell_0(\theta_0)}{\partial r_{i,j}^2} \right) = \frac{\partial^2}{\partial s^2} \left\{ \log \left( \tilde{\ell}_{i,j}(\gamma_0, s, Z_0) \right) \right\}_{s=r_{0,i,j}}, \quad 1 \le i < j \le k.$$

It is thus only necessary to check that for any $1 \le i < j \le k$,

$$s(i,j) = \mathbb{E} \left[ \frac{1}{\tilde{\ell}_{i,j}(\gamma_0, r_{0,i,j}, Z_0)^2} \left( \frac{\partial^2 \tilde{\ell}_{i,j}(\gamma_0, r_{0,i,j}, Z_0)}{\partial s^2} \times \tilde{\ell}_{i,j}(\gamma_0, r_{0,i,j}, Z_0) - \left( \frac{\partial \tilde{\ell}_{i,j}(\gamma_0, r_{0,i,j}, Z_0)}{\partial s} \right)^2 \right) \right] < 0.$$

Since $\tilde{\ell}_{i,j}(\gamma_0, r_{0,i,j}, Z_0)$ is the density of the pair $(Y_{i,0}, Y_{j,0})$ conditionally on $\mathcal{F}_{-1}$, evaluated at point $(Y_{i,0}, Y_{j,0})$, we have

$$\mathbb{E} \left[ \frac{\partial^2 \tilde{\ell}_{i,j}(\gamma_0, r_{0,i,j}, Z_0)}{\partial s^2} \right] = 0,$$

hence

$$s(i,j) = -\mathbb{E} \left[ \frac{1}{\tilde{\ell}_{i,j}(\gamma_0, r_{0,i,j}, Z_0)^2} \times \left( \frac{\partial \tilde{\ell}_{i,j}(\gamma_0, r_{0,i,j}, Z_0)}{\partial s} \right)^2 \right] \le 0.$$

However, if $s(i,j) = 0$, we have $\frac{\partial \tilde{\ell}_{i,j}(\gamma_0, r_{0,i,j}, Z_0)}{\partial s} = 0$ a.s. and when $Y_{i,0} = Y_{j,0} = 0$ we obtain

$$\frac{\partial}{\partial s} \left( s \mapsto \int_{-\infty}^{-\lambda_{i,0}} \Phi \left( -\frac{\lambda_{j,0} + sx}{\sqrt{1-s^2}} \right) \varphi(x) dx \right) \Big|_{s=r_{0,i,j}} = 0 \quad \text{a.s.}$$

But we already mention in Lemma 7.4 that such a derivative is always positive. Thus $s(i,j) < 0$ and $L_{2,2}$ is invertible. Thus, assumption **I7** is satisfied.
▷ We finally prove the following martingale difference property **I8**

$$ \forall t \in \mathbb{Z}, \mathbb{E}(\nabla_2 \ell_t(\theta_0) | \mathcal{F}_{t-1}) = 0. \quad (36) $$

Using the representation of $\ell_t(\theta_0)$ as a sum of terms depending on different parameters, we only have to prove that for $1 \le i < j \le k$

$$ \mathbb{E} \left( \frac{\partial \tilde{\ell}_{i,j}(\gamma_0, r_{0,i,j}, Z_t)}{\partial s} \mid \mathcal{F}_{t-1} \right) = 0. \quad (37) $$

Note that $\tilde{\ell}_{i,j}(\gamma_0, r_0, Z_t)$ is simply the conditional density of the pair $(Y_{i,t}, Y_{j,t})$ conditionally to $\mathcal{F}_{t-1}$. Equality (37) thus holds in a similar way as in the previous points, and assumption **I8** is satisfied.

One can finally state that $\sqrt{T}(\hat{\gamma} - \gamma_0, \hat{r} - r_0)$ has a Gaussian limiting distribution by plugging the result for estimator $\hat{\gamma}$ in equation (18). It yields, using the notation $m_\gamma(Z_t)$ defined previously

$$ \sqrt{T}(\hat{r} - r_0) = -\frac{L_{2,2}^{-1}}{\sqrt{T}} \sum_{t=p+1}^{T} \left\{ \nabla_2 \ell_t(\theta_0) - L_{i,2} J^{-1} \underbrace{\sum_{i=1}^{k} \hat{h}_{Y_{i,t}}(\lambda_{i,t}(\gamma_0)) \nabla_1 \lambda_{i,t}(\gamma_0)}_{\nabla_1 m_{\gamma_0}(Z_t)} \right\} + o_P(1). $$

Hence

$$ \begin{aligned} \sqrt{T}(\hat{\gamma} - \gamma_0, \hat{r} - r_0) &= \frac{1}{\sqrt{T}} \sum_{t=p+1}^{T} \left( J^{-1} \nabla_1 m_{\gamma_0}(Z_t), -L_{2,2}^{-1} (\nabla_2 \ell_t(\theta_0) - L_{1,2} J^{-1} \nabla_1 m_{\gamma_0}(Z_t)) \right) + o_P(1) \\ &= \left( \begin{array}{c|c} J^{-1} & 0 \\ \hline 0 & -L_{2,2}^{-1} \end{array} \right) \cdot \frac{1}{\sqrt{T}} \underbrace{\left( \begin{array}{c} \nabla_1 m_{\gamma_0}(Z_t) \\ \hline \nabla_2 \ell_t(\theta_0) - L_{1,2} J^{-1} \nabla_1 m_{\gamma_0}(Z_t) \end{array} \right)}_{W_t} + o_P(1). \end{aligned} $$

Yet, it has already been established that $(\nabla_1 m_{\gamma_0}(Z_t))_{t \in \mathbb{Z}}$ is a martingale difference series relatively to $(\mathcal{F}_t)_{t \in \mathbb{Z}}$. Then, by equality (36), $(W_t)_{t \in \mathbb{Z}}$ is also a martingale difference series relatively to $(\mathcal{F}_t)_{t \in \mathbb{Z}}$. Since we have already showed that $\mathbb{E}(\|W_0\|^2) < +\infty$, Theorem 7.13 then ensures the Gaussian limiting distribution of $\frac{W_t}{\sqrt{T}}$, and by Slutsky's lemma, the same conclusion holds for $\sqrt{T}(\hat{\gamma} - \gamma_0, \hat{r} - r_0)$.

**Lemma 7.7.** Take a and b two real numbers. Setting $\sigma = (1 + 4a^2)^{-1/2}$, we have the lower bound

$$ \int_{I_{Y_{i,0}-\lambda_{i,0}(\gamma)}} \Phi(ax+b)\varphi(x)dx \ge d^2\sigma \exp\left(-2b^2 - \frac{\lambda_{i,0}(\gamma)^2}{\sigma^2}\right) $$

where d is the positive constant mentioned in Lemma 7.8.

*Proof of Lemma 7.7.* By Lemma 7.8, there exists a positive constant $d > 0$ such that for all $x \in \mathbb{R}$

$$ \Phi(ax + b) \ge d \exp(-(ax + b)^2) \ge \exp(-2a^2x^2 - 2b^2). $$

We then get

$$ \begin{aligned} \int_{I_{Y_{i,0}-\lambda_{i,0}(\gamma)}} \Phi(ax+b)\varphi(x)dx &\ge d \exp(-2b^2) \int_{I_{Y_{i,0}-\lambda_{i,0}(\gamma)}} \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{x^2}{2\sigma^2}\right)dx \\ &\ge d \exp(-2b^2) \sigma \min\left\{\Phi\left(-\frac{\lambda_{i,0}(\gamma)}{\sigma}\right), 1-\Phi\left(-\frac{\lambda_{i,0}(\gamma)}{\sigma}\right)\right\} \\ &\ge d^2 \exp(-2b^2) \sigma \exp\left(-\frac{\lambda_{i,0}(\gamma)^2}{\sigma^2}\right). \end{aligned} $$

□
## 7.4 Proof of Proposition 4.1

Our aim here is to apply Theorem 7.14 from Giap and Van Quang (2016). For convenience, we also define $(Z_{j,t})_{t \in \mathbb{Z}}$ as independent copies of $(Z_{1,t})_{t \in \mathbb{Z}}$ for any negative $j$. We also set $Z := (Z_{j,t})_{(j,t) \in \mathbb{Z}^2}$.

The probability space $\Omega$ we work with here is $\mathbb{R}^{\mathbb{Z} \times \mathbb{Z}}$, embedded with the sigma field denoted $\mathbb{R}^{\mathbb{Z} \times \mathbb{Z}}$ and generated by the cylinders of the form

$$C = \prod_{j,t \in \mathbb{Z}} C_{j,t},$$

where $C_{j,t} \in \mathcal{B}(\mathbb{R})$, such that there exist $\mathcal{J}, \mathcal{T} \subset \mathbb{Z}$ finite and

$$\forall (j,t) \notin \mathcal{J} \times \mathcal{T}, C_{j,t} = \mathbb{R}.$$

The transformations of $\mathbb{R}^{\mathbb{Z} \times \mathbb{Z}}$ we consider here are

$$\tau_1: \mathbb{R}^{\mathbb{Z} \times \mathbb{Z}} \longrightarrow \mathbb{R}^{\mathbb{Z} \times \mathbb{Z}} \\ (x_{j,t})_{j,t \in \mathbb{Z}} \longmapsto (x_{j+1,t})_{j,t \in \mathbb{Z}}$$

$$\text{and} \quad \tau_2: \mathbb{R}^{\mathbb{Z} \times \mathbb{Z}} \longrightarrow \mathbb{R}^{\mathbb{Z} \times \mathbb{Z}} \\ (x_{j,t})_{j,t \in \mathbb{Z}} \longmapsto (x_{j,t+1})_{j,t \in \mathbb{Z}}.$$

Finally, the probability measure we consider is $\mathbb{P}_Z$

$$\forall A = (A_{j,t})_{j,t \in \mathbb{Z}} \in \mathbb{R}^{\mathbb{Z} \times \mathbb{Z}}, \mathbb{P}_Z(A) = \mathbb{P}\left(\bigcap_{j,t \in \mathbb{Z}} Z_{j,t}^{-1}(A_{j,t})\right).$$

▷ We first check that $\tau_1$ and $\tau_2$ preserve $\mathbb{P}_Z$. To this end, take any $A \in \mathbb{R}^{\mathbb{Z} \times \mathbb{Z}}$, we have

$$\mathbb{P}_Z(\tau_1^{-1}(A)) = \mathbb{P}\left(\bigcap_{j \in \mathbb{Z}} \bigcap_{t \in \mathbb{Z}} \left\{ Z_{j,t} \in A_{j+1,t} \right\} \right) \\ = \prod_{j \in \mathbb{Z}} \mathbb{P}(\bigcap_{t \in \mathbb{Z}} \left\{ Z_{j,t} \in A_{j+1,t} \right\})$$

by independence of the $(\mathcal{G}_j)_{j \in \mathbb{Z}}$

$$= \prod_{j \in \mathbb{Z}} \mathbb{P}(\bigcap_{t \in \mathbb{Z}} \{Z_{j+1,t} \in A_{j+1,t}\})$$

since the distribution of $(Z_{j,t})_{t \in \mathbb{Z}}$ does not depend on $j$

$$= \mathbb{P}_Z(A).$$

Hence $\tau_1$ preserves $\mathbb{P}_Z$. On the other hand

$$\mathbb{P}_Z(\tau_2^{-1}) = \mathbb{P}(\bigcap_{j \in \mathbb{Z}} \bigcap_{t \in \mathbb{Z}} \{Z_{j,t} \in A_{j,t+1}\}) \\ = \prod_{j \in \mathbb{Z}} \mathbb{P}(\bigcap_{t \in \mathbb{Z}} \{Z_{j,t} \in A_{j,t+1}\}) \\ = \prod_{j \in \mathbb{Z}} \mathbb{P}(\bigcap_{t \in \mathbb{Z}} \{Z_{j,t+1} \in A_{j,t+1}\})$$

by stationarity of $(Z_{j,t})_{t \in \mathbb{Z}}$

$$= \mathbb{P}_Z(A),$$

so $\tau_2$ also preserves $\mathbb{P}_Z$.

▷ Setting $h: (\mathbb{R}^d)^{\mathbb{Z} \times \mathbb{Z}} \to \mathbb{R}^d$ the mapping defined by $h(z) = z_{0,0}$, one can apply Theorem 7.14 to the mapping $f \circ h$ to get

$$\lim_{N,T \to +\infty} \frac{1}{NT} \sum_{j=1}^{N} \sum_{t=1}^{T} f \circ h \circ \tau_1^j \tau_2^t = \mathbb{E}_Z (f \circ h | \mathcal{I}_\tau) \quad \mathbb{P}_Z \text{ - a.s.}$$

where $\mathcal{I}_\tau = \mathcal{I}_{\tau_1} \cap \mathcal{I}_{\tau_2}$ and for $i \in \{1,2\}$, $\mathcal{I}_{\tau_i} = \{A \in \mathbb{R}^{\mathbb{Z} \times \mathbb{Z}} | \tau_i^{-1}(A) = A\}$.
▷ In order to obtain the desired result, it is then sufficient to prove that $\tau_1$ is ergodic, i.e.

$$ \forall A \in \mathcal{I}_{\tau_1}, \mathbb{P}_Z(A) \in \{0, 1\}. $$

We actually show that for any $A \in \mathcal{I}_{\tau_1}$

$$ \lim_{n \to +\infty} \mathbb{P}_Z(A \cap \tau_1^{-n}(A)) = \mathbb{P}_Z(A)^2. \quad (39) $$

First, this latter equality is satisfied for any cylinder of the form

$$ B = \prod_{j,t \in \mathbb{Z}} B_{j,t}, $$

with for all $|j| > J$

$$ B_{j,t} = \mathbb{R}, \forall t \in \mathbb{Z}, $$

for some fixed $J \in \mathbb{N}$. Indeed, for all $n \ge 2J + 1$, we have

$$ \begin{aligned} \mathbb{P}_Z(B \cap \tau_1^{-n}(B)) &= \mathbb{P} \left( \bigcap_{j \in \mathbb{Z}} \cap_{t \in \mathbb{Z}} \{ Z_{j,t} \in B_{j,t} \} \cap \bigcap_{j \in \mathbb{Z}} \cap_{t \in \mathbb{Z}} \{ Z_{j+n,t} \in B_{j,t} \} \right) \\ &= \mathbb{P} \left( \bigcap_{j=-J}^{J} \underbrace{\cap_{t \in \mathbb{Z}} \{ Z_{j,t} \in B_{j,t} \}}_{\in \mathcal{G}_j} \cap \bigcap_{j'=n-J}^{n+J} \underbrace{\cap_{t \in \mathbb{Z}} \{ Z_{j',t} \in B_{j'-n,t} \}}_{\in \mathcal{G}_{j'}} \right). \end{aligned} $$

Yet, the indices $j$ and $j'$ are necessarily different since $n \ge 2J + 1$, so by independence of the $(\mathcal{G}_j)_{j \in \mathbb{Z}}$

$$ \begin{aligned} \mathbb{P}_Z(B \cap \tau_1^{-n}(B)) &= \mathbb{P}_Z(B) \times \prod_{j'=n-J}^{n+J} \mathbb{P} \left( \cap_{t \in \mathbb{Z}} \{ Z_{j',t} \in B_{j'-n,t} \} \right) \\ &= \mathbb{P}_Z(B) \times \prod_{j=-J}^{J} \mathbb{P} \left( \cap_{t \in \mathbb{Z}} \{ Z_{j+n,t} \in B_{j,t} \} \right) \\ &= \mathbb{P}_Z(B) \times \prod_{j=-J}^{J} \mathbb{P} \left( \cap_{t \in \mathbb{Z}} \{ Z_{j,t} \in B_{j,t} \} \right) \end{aligned} $$

since the distribution of $(Z_{j,t})_{t \in \mathbb{Z}}$ does not depend on $j$

$$ = \mathbb{P}_Z(B)^2. \quad (40) $$

We now prove equality (39), and take any $A \in \mathcal{I}_{\tau_1}$ and $\varepsilon > 0$. There exists a cylinder $B$ such that $\mathbb{P}_Z(A\Delta B) < \varepsilon$ and we have

$$ \begin{aligned} |\mathbb{P}_Z(A \cap \tau_1^{-n}(A)) - \mathbb{P}_Z(A)^2| &\le |\mathbb{P}_Z(A \cap \tau_1^{-n}(A)) - \mathbb{P}_Z(B \cap \tau_1^{-n}(A))| \\ &\quad + |\mathbb{P}_Z(B \cap \tau_1^{-n}(A)) - \mathbb{P}_Z(B \cap \tau_1^{-n}(B))| \\ &\quad + |\mathbb{P}_Z(B \cap \tau_1^{-n}(B)) - \mathbb{P}_Z(B)^2| \\ &\quad + |\mathbb{P}_Z(B)^2 - \mathbb{P}_Z(A)^2|. \end{aligned} $$

However

$$ |\mathbb{P}_Z(B)^2 - \mathbb{P}_Z(A)^2| \le 2 \times |\mathbb{P}_Z(B) - \mathbb{P}_Z(A)| \le 2\varepsilon $$

and for $n$ large enough, we have according to equality (40)

$$ |\mathbb{P}_Z(B \cap \tau_1^{-n}(B)) - \mathbb{P}_Z(B)^2| = 0. $$

For such a $n$, we also have since $A \in \mathcal{I}_{\tau_1}$

$$ \begin{aligned} |\mathbb{P}_Z(A \cap \tau_1^{-n}(A)) - \mathbb{P}_Z(B \cap \tau_1^{-n}(A))| &= |\mathbb{P}_Z(A) - \mathbb{P}_Z(B \cap A)| \\ &= \mathbb{P}_Z(A \setminus B) \\ &< \varepsilon. \end{aligned} $$
Finally, we have

$$
\begin{align*}
|\mathbb{P}_Z (B \cap \tau_1^{-n}(A)) - \mathbb{P}_Z (B \cap \tau_1^{-n}(B))| &\le \mathbb{P}_Z (\tau_1^{-n}(A)\Delta\tau_1^{-n}(B)) \\
&= \mathbb{P}_Z (\tau_1^{-n} (A\Delta B)) \\
&= \mathbb{P}_Z (A\Delta B) \quad \text{since } \tau_1 \text{ preserves } \mathbb{P}_Z \\
&< \varepsilon.
\end{align*}
$$

Thus, for $n \in \mathbb{N}$ large enough

$$
|\mathbb{P}_Z (A \cap \tau_1^{-n}(A)) - \mathbb{P}_Z (A)^2| < 4\varepsilon,
$$

hence

$$
\lim_{n \to +\infty} \mathbb{P}_Z (A \cap \tau_1^{-n}(A)) = \mathbb{P}_Z (A) = \mathbb{P}_Z (A)^2.
$$

We then easily deduce that $\mathbb{P}_Z(A) \in \{0,1\}$, and we have the desired result.

7.5 Proof of Proposition 4.2

As in the proof of Proposition 4.1, we may assume that $Z_{j,t}$ is also defined for a non-positive integer $j$.
To this end, it is simply necessary to define $(V_j, \kappa_{j,t}, \varepsilon_{j,t})_{t \in \mathbb{Z}}$ for non-positive indices $j$, independently
across the $j$ and with the same distribution as $(V_1, \kappa_{1,t}, \varepsilon_{1,t})_{t \in \mathbb{Z}}$. Set $\tilde{\kappa}_{j,t} = (\eta_t, \kappa_{j,t}, \varepsilon_{j,t})$ for $(j,t) \in \mathbb{N} \times \mathbb{Z}$.
It is quite clear that there exists a measurable mapping $\bar{G} : \mathbb{R}^{d_2} \times (\mathbb{R}^{d_4+d_6+k})^{\mathbb{N}} \to \mathbb{R}$ such that

$$
Z_{j,t} = \bar{G} \left( V_j, (\tilde{\kappa}_{j,t-s})_{s \ge 0} \right).
$$

From our independence conditions, it is straightforward to show that the probability distribution of
$Z := (Z_{j,t})_{(j,t) \in \mathbb{Z}^2}$ is invariant under the two shifts $\tau_1$ and $\tau_2$ defined in the proof of Proposition 4.1.
We are going to show that for all Borel subset $A$ of $\mathbb{R}^{\mathbb{Z} \times \mathbb{Z}}$, then

$$
\limsup_{m \to \infty} \limsup_{n \to \infty} \mathbb{P}_Z (A \cap \tau_1^{-m} \tau_2^{-n} A) = \mathbb{P}_Z (A)^2. \quad (41)
$$

Applying (41) to $A \in \mathcal{I}_{\tau_1} \cap \mathcal{I}_{\tau_2}$ an invariant set, we will get $\mathbb{P}_Z(A) \in \{0,1\}$. We will prove (41) when $A$
is cylinder set. The rest of the proof is similar to that of Proposition 4.1. For $(j,t) \in \mathbb{Z}^2$ and $v_j \in \mathbb{R}^{d_2}$,
set

$$
Z_{j,t}(v_j) = \bar{G} \left( v_j, (\tilde{\kappa}_{j,t-s})_{s \ge 0} \right).
$$

If $A = \prod_{(j,t) \in \mathbb{Z}^2} A_{j,t}$ with $A_{j,t} = \mathbb{R}$ except if $(j,t) \in J \times K$, for some $J$ and $K$ some finite subsets of $\mathbb{Z}$, we have

$$
\mathbb{P}_Z (A \cap \tau_1^{-m} \tau_2^{-n} A) = \mathbb{P} (\cap_{(j,t) \in J \times K} \{Z_{j,t} \in A_{j,t}, Z_{j+m, t+n} \in A_{j,t}\}).
$$

Conditioning with respect to $V_j = v_j$ for $j \in \mathbb{Z}$, we get for a fixed value of $m$,

$$
\begin{align*}
& \lim_{n \to \infty} \mathbb{P} \left( \cap_{(j,t) \in J \times K} \{ Z_{j,t}(v_j) \in A_{j,t}, Z_{j+m,t+n}(v_{j+m}) \in A_{j,t} \} \right) \\
& = \mathbb{P} \left( \cap_{(j,t) \in J \times K} \{ Z_{j,t}(v_j) \in A_{j,t} \} \right) \mathbb{P} \left( \cap_{(j,t) \in J \times K} \{ Z_{j,t}(v_{j+m}) \in A_{j,t} \} \right).
\end{align*}
$$

This property is due to the mixing property of the Bernoulli shifts. See for instance Samorodnitsky
et al. (2016), Corollary 2.2.5. Indeed, setting $J_m = J \cup \{m + j : j \in J\}$, $x_{J_m} = (x_j)_{j \in J_m}$ for $x \in \mathbb{R}^Z$,
and for $t \in \mathbb{Z}$, $Z_{J_m, t} = (Z_{j,t}(v_j))_{j \in J_m}$, $B_{j,t} = A_{j-m, t}$ if $j-m \in J$ and $B_{j,t} = \mathbb{R}$ if $j-m \notin J$, we have
the equality

$$
\cap_{(j,t) \in J \times K} \{ Z_{j,t}(v_j) \in A_{j,t}, Z_{j+m,t+n}(v_{j+m}) \in A_{j,t} \} = \cap_{t \in K} \left\{ Z_{J_m, t} \in \prod_{j \in J_m} A_{j,t}, Z_{J_m, t+n} \in \prod_{j \in J_m} B_{j,t} \right\}.
$$

But since $Z_{J_m, t} \in \sigma(\tilde{\kappa}_{J_m, t-s} : s \ge 0)$, one can use the mixing property for this multivariate Bernoulli shift.
Now if $m$ is large enough, we have $j \neq j' + m$ for any pair $(j, j') \in J^2$. One can integrate the previous equality with respect to the distribution of the $V_j$s and then use the dominated convergence theorem as well as the independence properties to get

$$ \lim_{n \to \infty} \mathbb{P} \left( \bigcap_{(j,t) \in J \times K} \{ Z_{j,t} \in A_{j,t}, Z_{j+m,t+n} \in A_{j,t} \} \right) \\ = \mathbb{P} \left( \bigcap_{(j,t) \in J \times K} \{ Z_{j,t} \in A_{j,t} \} \right)^2. $$

This yields to (41) and completes the proof. $\square$

## 7.6 Proof of Theorem 4.3

First we note that under Assumption C1 or C'1 and Assumption C2, Corollary 3.1 guarantees the existence of a unique stationary solution $(Y_{j,t})_{t \in \mathbb{Z}}$ to the recursions.

All the asymptotic results follow by checking the assumptions of Theorem 7.15 and Theorem 7.17 and their corollaries. Similar assumptions have been already checked for proving Theorem 3.2 and no new arguments are necessary. Details are then omitted. $\square$

## 7.7 Additional asymptotic results for general M-estimators

This section consists in a collection of various results that are used in section 6. Proofs are omitted here.

We begin with a useful inequality about the CDF $\Phi$ of a Gaussian distribution $N(0, 1)$

**Lemma 7.8.** There exist two positive real numbers $c$ and $d$ such that for all $x \in \mathbb{R}$

$$ \min \{\Phi(x), 1 - \Phi(x)\} \geq \frac{c}{1 + |x|} e^{-x^2/2} \geq d e^{-x^2}. $$

The following result is derived from a more general result in Pfanzagl (1969), and is about consistency. We propose an alternative proof for our setup.

**Theorem 7.9.** Let $(Z_t)_{t \in \mathbb{Z}}$ be a stochastic process valued in a measurable space $\mathcal{Z}$, whose distribution depends on an unknown parameter $\theta_0 \in \Theta$, where $\Theta$ is a set of $\mathbb{R}^n$, $n \in \mathbb{N}^*$.

For all $\theta \in \Theta$, we consider a measurable application $m_\theta: \mathcal{Z} \to \mathbb{R}$.

Let us define the assumptions:

**H1** $(Z_t)_{t \in \mathbb{Z}}$ is ergodic.

**H2** $\Theta$ is a compact set.

**H3** For all $z \in \mathcal{Z}$, the application $\theta \mapsto m_\theta(z)$ is continuous and $\mathbb{E} \left( \sup_{\theta \in \Theta} |m_\theta(Z_0)| \right) < +\infty$.

**H4** For all $\theta \in \Theta$, $\mathbb{E}(m_\theta(Z_0)) \leq \mathbb{E}(m_{\theta_0}(Z_0))$ with equality if and only if $\theta = \theta_0$.

If assumptions **H1** to **H4** are satisfied, then

$$ \hat{\theta}_T := \underset{\theta \in \Theta}{\operatorname{argmax}} \sum_{t=1}^{T} m_{\theta}(Z_t) \underset{T \to +\infty}{\longrightarrow} \theta_0 \quad \text{a.s.} $$

The next result deals with asymptotic normality, and is very similar to Theorem 5.6.1 in Straumann (2006). Here, $||.||$ denotes any norm on the considered space, for example the infinite norm. Furthermore, for a generic function $f: \mathbb{R}^n \to \mathbb{R}$, we denote $\nabla f$ the gradient vector of $f$ and $\nabla^2 f$ its hessian matrix.
**Theorem 7.10.** Let us consider a process $(Z_t)_{t \in \mathbb{Z}}$ satisfying conditions **H1** to **H4** mentioned in theorem 7.9, and consider $(\mathcal{F}_t)_{t \in \mathbb{Z}}$ a filtration. We assume that for all $z \in \mathcal{Z}$, the application $\theta \mapsto m_\theta(z)$ is of class $\mathcal{C}^2$ on $\Theta$.

Assume also that

**H5** $\theta_0$ belongs to the interior of $\Theta$.

**H6** $\mathbb{E} (\|\nabla m_{\theta_0}(Z_0)\|^2) < +\infty$.

**H7** $\mathbb{E} \left( \sup_{\theta \in \Theta} \|\nabla^2 m_\theta(Z_0)\| \right) < +\infty$.

**H8** The matrix $\mathbb{E} (\nabla^2 m_{\theta_0}(Z_0))$ is invertible.

**H9** For all $t \in \mathbb{Z}$, $\mathbb{E}(\nabla m_{\theta_0}(Z_t) | \mathcal{F}_{t-1}) = 0$, i.e. $(\nabla m_{\theta_0}(Z_t))_{t \in \mathbb{Z}}$ is a martingale difference series relatively to $(\mathcal{F}_{t-1})_{t \in \mathbb{Z}}$.

Then, $\hat{\theta}_T := \underset{\theta \in \Theta}{\operatorname{argmax}} \sum_{t=1}^T m_\theta(Z_t) \underset{n \to +\infty}{\longrightarrow} \theta_0$ almost surely, and we have the stochastic expansion

$$ \mathbb{E} (\nabla^2 m_{\theta_0}(Z_0)) \cdot \sqrt{T} (\hat{\theta}_T - \theta_0) = - \frac{1}{\sqrt{T}} \sum_{t=1}^T \nabla m_{\theta_0}(Z_t) + o_P(1). $$

This leads to the following asymptotic normality

$$ \sqrt{n} (\hat{\theta}_n - \theta_0) \stackrel{\mathcal{L}}{\longrightarrow} \mathcal{N} \left( 0, \mathbb{E} \left( \nabla^2 m_{\theta_0}(Z_0) \right)^{-1} \cdot \operatorname{Var} \left( \nabla m_{\theta_0}(Z_0) \right) \cdot \mathbb{E} \left( \nabla^2 m_{\theta_0}(Z_0) \right)^{-1} \right). $$

Proposition 8 of Z.-M. Debaly and Truquet (2023) states the following refinement of Theorem 7.9

**Proposition 7.11.** Take $\Gamma$ and $\mathcal{R}$ two compact sets of $\mathbb{R}^p$ and $\mathbb{R}^q$ respectively, and denote $\Theta = \Gamma \times \mathcal{R}$.

Let $(Z_t)_{t \in \mathbb{Z}}$ be a stochastic process valued in a measurable space $\mathcal{Z}$, whose distribution depends on an unknown parameter $\theta_0 = (\gamma_0, r_0) \in \Theta$.

For all $\theta = (\gamma, r) \in \Theta$, we consider a measurable application $m_\theta : \mathcal{Z} \longrightarrow \mathbb{R}$.

Let us define the assumptions:

**I1** $(Z_t)_{t \in \mathbb{Z}}$ is ergodic.

**I2** For all $z \in \mathcal{Z}$, the application $\theta \mapsto m_\theta(z)$ is continuous and $\mathbb{E} \left( \sup_{\theta \in \Theta} |m_\theta(Z_0)| \right) < +\infty$.

**I3** For all $r \in \mathcal{R}$, $\mathbb{E} (m_{(\gamma_0, r)}(Z_0)) \le \mathbb{E} (m_{(\gamma_0, r_0)}(Z_0))$ with equality if and only if $r = r_0$.

Finally, we suppose that for any $T \in \mathbb{N}^*$, we have a strongly consistent estimator $\hat{\gamma}_T$ of $\gamma_0$

$$ \hat{\gamma}_T \underset{T \to +\infty}{\longrightarrow} \gamma_0 \quad \text{a.s.} $$

Then, if assumptions **I1** to **I3** are satisfied, we have

$$ \hat{r}_T := \underset{r \in \mathcal{R}}{\operatorname{argmax}} \sum_{t=1}^T m_{(\hat{\gamma}_T, r)}(Z_t) \underset{T \to +\infty}{\longrightarrow} r_0 \quad \text{a.s.} $$

It is also possible to obtain a refinement of Theorem 7.10, that we detail in the following result.

Here, for a generic function

$$ f: \mathbb{R}^p \times \mathbb{R}^q \longrightarrow \mathbb{R} \\ (x, y) \longmapsto f(x, y) $$

we denote by $\nabla_1 f$ and $\nabla_2 f$ the gradients of function $f$ obtained when we differentiate respectively with respect to $x$ and $y$. Similarly, we denote $\nabla_1^2 f$ and $\nabla_2^2 f$ the corresponding Hessian matrices, and we also denote

$$ \nabla_{1,2} f = \left( \frac{\partial^2 f}{\partial x_j \partial y_i} \right)_{1 \le i \le q, 1 \le j \le p}. $$
**Proposition 7.12.** Take $\Gamma$ and $\mathcal{R}$ two compact sets of $\mathbb{R}^p$ and $\mathbb{R}^q$ respectively, and denote $\Theta = \Gamma \times \mathcal{R}$.
Let $(Z_t)_{t \in \mathbb{Z}}$ be a stochastic process satisfying conditions **I1** to **I3** mentioned in Proposition 7.11, and
consider $(\mathcal{F}_t)_{t \in \mathbb{Z}}$ a filtration. Assume that for all $z \in \mathbb{Z}$, the application $\theta \mapsto m_\theta(z)$ is of class $C^2$ on
$\Theta$. Assume also that

**I4** $\theta_0$ belongs to the interior of $\Theta$.

**I5** $\mathbb{E} (\|\nabla_2 m_{\theta_0}(Z_0)\|^2) < +\infty$.

**I6** $\mathbb{E} \left( \sup_{\theta \in \Theta} \| \nabla_{1,2} m_{\theta}(Z_0) \| \right) < +\infty$ and $\mathbb{E} \left( \sup_{\theta \in \Theta} \| \nabla_2^2 m_{\theta}(Z_0) \| \right) < +\infty$.

**I7** The matrix $\mathbb{E} (\nabla_2^2 m_{\theta_0}(Z_0))$ is invertible.

**I8** For all $t \in \mathbb{Z}$, $\mathbb{E}(\nabla_2 m_{\theta_0}(Z_t) | \mathcal{F}_{t-1}) = 0$, i.e. $(\nabla_2 m_{\theta_0}(Z_t))_{t \in \mathbb{Z}}$ is a martingale difference series relatively to $(\mathcal{F}_{t-1})_{t \in \mathbb{Z}}$.

Finally, we suppose that for any $T \in \mathbb{N}^*$, we have a strongly consistent estimator $\hat{\gamma}_T$ of $\gamma_0$ such that

$$ \sqrt{T}(\hat{\gamma}_T - \gamma_0) = O_P(1). $$

Then $\hat{r}_T := \underset{r \in \mathcal{R}}{\operatorname{argmax}} \sum_{t=1}^T m_{(\hat{\gamma}_T, r)}(Z_t) \underset{T \to +\infty}{\longrightarrow} r_0$ almost surely, and we have the stochastic expansion

$$ \mathbb{E} (\nabla_2^2 m_{\theta_0}(Z_0)) \cdot \sqrt{T} (\hat{r}_T - r_0) = - \frac{1}{\sqrt{T}} \sum_{t=1}^T \nabla_2 m_{\theta_0}(Z_t) - \mathbb{E} (\nabla_{1,2} m_{\theta_0}) \sqrt{T} (\hat{\gamma}_T - \gamma_0) + o_P(1). $$

The previous theorems dealing about asymptotic normality are largely based on the well-known Ibragimov-Billingsley's theorem for martingale differencies below. See Billingsley (2013), Theorem 35.12, for a proof. This theorem extends naturally to multivariate random variables.

**Theorem 7.13.** Consider an array $(Y_{k,t})_{k,t \in \mathbb{N}^*}$ an array of real random variables such that for each $k \in \mathbb{N}^*$, $(Y_{k,t})_{t \in \mathbb{N}^*}$ is a difference martingale series with respect to a filtration $(\mathcal{G}_{k,t})_{t \in \mathbb{N}}$ (with $\mathcal{G}_{k,0} = \{\emptyset, \Omega\}$).

Suppose the $Y_{k,t}$'s have second moments, and set

$$ \sigma_{k,t}^2 = \mathbb{E} (Y_{k,t}^2 | \mathcal{G}_{k,t-1}). $$

Assume that

$$ \sum_{t=1}^{+\infty} \sigma_{k,t}^2 \xrightarrow{\mathbb{P}} \sigma^2, $$

where $\sigma$ is a positive constant, and that

$$ \sum_{t=1}^{+\infty} \mathbb{E} \left( Y_{k,t}^2 \mathbb{1}_{|Y_{k,t}| \ge \varepsilon} \right) \underset{k \to +\infty}{\longrightarrow} 0 $$

for each $\varepsilon > 0$. Then

$$ \sum_{t=1}^{+\infty} Y_{k,t} \xrightarrow{\mathcal{L}} \mathcal{N}(0, \sigma^2) $$

when $k \to +\infty$.

Let us finally recall Lemma 3.1 from Giap and Van Quang (2016), which is a fundamental result in our paper.
**Theorem 7.14.** Consider $(\mathcal{X}, \|\cdot\|)$ a Banach space and let $k$ be a positive integer, and $\theta_1, \dots, \theta_k$ be commutative measure preserving transformations of a probability space $(\Omega, \mathcal{A}, \mathbb{P})$.

Then for any measurable application $f: \Omega \to \mathcal{X}$ such that $\mathbb{E}(\|f\|(\log^+(\|f\|)^{k-1}) < +\infty$, the multiple sequence

$$A_{m_1, \dots, m_k} f = \frac{1}{m_1 \dots m_k} \sum_{i_1=0}^{m_1-1} \dots \sum_{i_k=0}^{m_k-1} f \left( \theta_1^{i_1} \dots \theta_k^{i_k} \right)$$

converges almost surely to $\mathbb{E}(f | \mathcal{I}_{\theta})$ as $\min(m_1, \dots, m_k) \to +\infty$, where $\mathcal{I}_{\theta} = \cap_{i=1}^k \mathcal{I}_{\theta_i}$ and $\mathcal{I}_{\theta_i} = \{A \in \mathcal{A} : \theta_i^{-1}(A) = A\}$.

Moreover, if $\theta_s$ is ergodic for some $s$ belonging to $\{1, \dots, k\}$, $\mathbb{E}(f | \mathcal{I}_{\theta}) = \mathbb{E}(f)$ almost surely.

## 7.8 Asymptotic results for M-estimators with panel data

We now state several general results relative to M-estimators for panel data. The two first ones deal with the consistency, while the others allow us to obtain central limit theorems.

**Theorem 7.15.** Let $(Z_{j,t})_{j \in \mathbb{N}^*, t \in \mathbb{Z}}$ be a family of random variables valued in $\mathcal{Z} = \mathbb{R}^{d'}$ for some positive integer $d'$, satisfying (11) and either C1 or C'1.

Additionally, for a parameter $\theta$ belonging to a compact set $\Theta \subset \mathbb{R}^{d_{\Theta}}$, let us consider a measurable function $m_{\theta}: \mathcal{Z} \to \mathbb{R}$, and assume that:

**P1** For all $z \in \mathcal{Z}$, the application $\theta \mapsto m_{\theta}(z)$ is continuous.

**P2** There exists $\delta > 0$ such that $\mathbb{E} \left( \sup_{\theta \in \Theta} |m_{\theta}(Z_{1,1})|^{1+\delta} \right) < +\infty$.

**P3** There exists $\theta_0 \in \Theta$ such that for all $\theta \in \Theta$, $\mathbb{E}(m_{\theta}(Z_{1,1})) \le \mathbb{E}(m_{\theta_0}(Z_{1,1}))$, with equality if and only if $\theta = \theta_0$.

Then

$$\hat{\theta}_{n,T} := \underset{\theta \in \Theta}{\operatorname{argmax}} \sum_{t=1}^{T} \sum_{j=1}^{n} m_{\theta}(Z_{j,t}) \underset{n,T \to +\infty}{\longrightarrow} \theta_0 \quad \text{a.s.}$$

It is possible to refine Theorem 7.15 in the following case.

**Corollary 7.16.** Let $(Z_{j,t})_{j \in \mathbb{N}^*, t \in \mathbb{Z}}$ be a family of random variables valued in $\mathcal{Z} = \mathbb{R}^{d'}$ for some positive integer $d'$, satisfying (11) and either C1 or C'1. We consider two compact sets $\Gamma$ and $\mathcal{R}$ of $\mathbb{R}^{d_{\Gamma}}$ and $\mathbb{R}^{d_{\mathcal{R}}}$ respectively, and we denote $\Theta = \Gamma \times \mathcal{R}$.

Additionally, for a parameter $\theta = (\gamma, r) \in \Theta$, we consider a measurable function $m_{\theta}: \mathcal{Z} \to \mathbb{R}$, and assume that:

**Q1** For all $z \in \mathcal{Z}$, the application $\theta \mapsto m_{\theta}(z)$ is continuous.

**Q2** $\mathbb{E} \left( \sup_{\theta \in \Theta} |m_{\theta}(Z_{1,1})|^{1+\delta} \right) < +\infty$ for some $\delta > 0$.

**Q3** There exists $\theta_0 = (\gamma_0, r_0) \in \Theta$ such that for all $r \in \mathcal{R}$, $\mathbb{E}(m_{\gamma_0, r}(Z_{1,1})) \le \mathbb{E}(m_{\gamma_0, r_0}(Z_{1,1}))$, with equality if and only if $r = r_0$.

If $\hat{\gamma}_{n,T}$ is a statistic such that $\hat{\gamma}_{n,T} \underset{n,T \to +\infty}{\longrightarrow} \gamma_0$ almost surely, then

$$\hat{r}_{n,T} := \underset{r \in \mathcal{R}}{\operatorname{argmax}} \sum_{t=1}^{T} \sum_{j=1}^{n} m_{\hat{\gamma}_{n,T}}(Z_{j,t}) \underset{n,T \to +\infty}{\longrightarrow} r_0 \quad \text{a.s.}$$
In what follows, for any function $f_\theta$ valued in $\mathbb{R}$, depending on a parameter $\theta \in \mathbb{R}^{d_\Theta}$, we denote respectively $\nabla f_\theta$ and $\nabla^2 f_\theta$ the gradient and hessian matrix of $f_\theta$ when we differentiate with respect to $\theta$. We also denote by $\mathcal{A} \lor \mathcal{B}$ the sigma-field generated by two sigma-fields $\mathcal{A}$ and $\mathcal{B}$.

Furthermore, $\| \cdot \|$ denotes any norm on the considered vector space.

**Theorem 7.17.** Let $(Z_{j,t})_{j \in \mathbb{N}^*, t \in \mathbb{Z}}$ be a family of random variables valued in $\mathcal{Z} = \mathbb{R}^{d'}$ for some positive integer $d'$, satisfying (11) and either **C1** or **C'1**.

Additionally, for a parameter $\theta$ belonging to a compact set $\Theta \subset \mathbb{R}^{d_\Theta}$, let us consider a measurable function $m_\theta: \mathcal{Z} \to \mathbb{R}$ such that assumptions **P1** to **P3** are satisfied.

Assume furthermore that

**P4** $\theta_0$ belongs to the interior of $\Theta$.

**P5** For all $z \in \mathcal{Z}$, $\theta \mapsto m_\theta(z)$ is of class $C^2$.

**P6** $\mathbb{E} (\|\nabla m_{\theta_0}(Z_{1,1})\|^{2+\delta}) < +\infty$ for some $\delta > 0$.

**P7** $\mathbb{E} \left( \sup_{\theta \in \Theta} \|\nabla^2 m_\theta(Z_{1,1})\|^{1+\delta} \right) < +\infty$ for some $\delta > 0$.

**P8** $\mathbb{E} (\nabla^2 m_{\theta_0}(Z_{1,1}))$ is invertible.

**P9** Denoting for $t \in \mathbb{Z}$, $\mathcal{H}_t = \cup_{\ell \ge 1} \mathcal{F}_{\ell,t}$ where the $\mathcal{F}_{\ell,t}$s are defined in (12), then $\mathbb{E} (\nabla m_{\theta_0}(Z_{j,t}) | \mathcal{H}_{t-1}) = 0$ for $t \in \mathbb{Z}$ and $j \ge 1$.

Then $\hat{\theta}_{n,T} := \underset{\theta \in \Theta}{\operatorname{argmax}} \sum_{t=1}^T \sum_{j=1}^n m_\theta(Z_{j,t}) \underset{n,T \to +\infty}{\longrightarrow} \theta_0$ almost surely and we have the stochastic expansion

$$ \sqrt{nT} \mathbb{E} (\nabla^2 m_{\theta_0}(Z_{1,1})) \cdot (\hat{\theta}_{n,t} - \theta_0) = - \frac{1}{\sqrt{nT}} \sum_{t=1}^T \sum_{j=1}^n \nabla m_{\theta_0}(Z_{j,t}) + o_P(1). $$

when $n, T \to \infty$. Therefore, we have the convergence in distribution

$$ \sqrt{nT} (\hat{\theta}_{n,T} - \theta_0) \xrightarrow{\mathcal{L}} \mathcal{N} \left( 0, \mathbb{E} \left( \nabla^2 m_{\theta_0}(Z_{1,1}) \right)^{-1} \cdot \mathrm{Var} \left( \nabla m_{\theta_0}(Z_{1,1}) \right) \mathbb{E} \left( \nabla^2 m_{\theta_0}(Z_{1,1}) \right)^{-1} \right) $$

when $n, T \to \infty$.

Similarly to Corollary 7.16, it is possible to refine Theorem 7.17 in the following case.

Just like in section 3, for a function $f_{(\gamma,r)}$ valued in $\mathbb{R}$, depending on a parameter $(\gamma, r) \in \Gamma \times \mathcal{R} \subset \mathbb{R}^{d_\Gamma} \times \mathbb{R}^{d_\Gamma}$, we denote $\nabla_1 f$ and $\nabla_2 f$ the gradients obtained by differentiating $f$ with respect to $\gamma$ or $r$ respectively. We denote in a similar way $\nabla_1^2$ and $\nabla_2^2$ the corresponding hessian matrices and

$$ \nabla_{1,2} f = \left( \frac{\partial^2 f}{\partial \gamma \partial r_i} \right)_{1 \le i \le d_{\mathcal{R}}, 1 \le j \le d_{\Gamma}} \in \mathbb{R}^{d_{\mathcal{R}} \times d_{\Gamma}}. $$

**Corollary 7.18.** Let $(Z_{j,t})_{j \in \mathbb{N}^*, t \in \mathbb{Z}}$ be a family of random variables valued in $\mathcal{Z} = \mathbb{R}^{d'}$ for some positive integer $d'$, satisfying (11) and either **C1** or **C'1**. We consider two compact sets $\Gamma$ and $\mathcal{R}$ of $\mathbb{R}^{d_\Gamma}$ and $\mathbb{R}^{d_\Gamma}$ respectively, and we denote $\Theta = \Gamma \times \mathcal{R}$.

Additionally, for a parameter $\theta = (\gamma, r) \in \Theta$, we consider a measurable function $m_\theta: \mathcal{Z} \to \mathbb{R}$ such that assumptions **Q1** to **Q3** are satisfied.

Assume furthermore that

**Q4** $\theta_0$ belongs to the interior of $\Theta$.

**Q5** For all $z \in \mathcal{Z}$, $\theta \mapsto m_\theta(Z_{1,1})$ is of class $C^2$.

**Q6** $\mathbb{E} (\|\nabla_2 m_{\theta_0}(Z_{1,1})\|^{2+\delta}) < +\infty$ for some $\delta > 0$.
$$ \mathbf{Q7} \quad \mathbb{E} \left( \sup_{\theta \in \Theta} \| \nabla_{1,2} m_{\theta}(Z_{1,1}) \|^{1+\delta} \right) < +\infty \text{ for some } \delta > 0. $$

$$ \mathbf{Q8} \quad \mathbb{E} \left( \sup_{\theta \in \Theta} \| \nabla_2^2 m_{\theta}(Z_{1,1}) \|^{1+\delta} \right) < +\infty \text{ for some } \delta > 0. $$

$$ \mathbf{Q9} \quad \mathbb{E} (\nabla_2^2 m_{\theta_0}(Z_{1,1})) \text{ is invertible.} $$

**Q10** Denoting for $t \in \mathbb{Z}$, $\mathcal{H}_t = \bigvee_{\ell \ge 1} \mathcal{F}_{\ell,t}$ where the $\mathcal{F}_{\ell,t}$s are defined in (12), then $\mathbb{E}(\nabla_2 m_{\theta_0}(Z_{j,t}) | \mathcal{H}_{t-1}) = 0$ for $t \in \mathbb{Z}$ and $j \ge 1$.

Finally, we suppose that for any $n, T$ in $\mathbb{N}^*$, we have a strongly consistent estimator $\hat{\gamma}_{n,T}$ of $\gamma_0$ such that

$$ \sqrt{nT} (\hat{\gamma}_{n,T} - \gamma_0) = o_{\mathbb{P}}(1). $$

when $n, T \to \infty$.

Then $\hat{r}_{n,T} := \underset{r \in \mathcal{R}}{\operatorname{argmax}} \sum_{t=1}^T \sum_{j=1}^n m_{(\hat{\gamma}_{n,T},r)}(Z_{j,t}) \underset{n,T \to +\infty}{\longrightarrow} r_0$ almost surely and we have the stochastic expansion

$$ \sqrt{nT} \mathbb{E} (\nabla_2^2 m_{\theta_0}(Z_{1,1})) (\hat{r}_{n,T} - r_0) = - \frac{1}{\sqrt{nT}} \sum_{t=1}^T \sum_{j=1}^n \nabla_2 m_{\theta_0}(Z_{j,t}) - \mathbb{E} (\nabla_{1,2} m_{\theta_0}(Z_{1,1})) \sqrt{nT} (\hat{\gamma}_{n,T} - \gamma_0) + o_{\mathbb{P}}(1). $$

**Note.** Once again, assumption **Q8** implies assumption **Q1**.

## 7.9 Proof of Theorem 7.15

The proof of Theorem 7.15 is actually an adaptation of the proof of Theorem 7.9 and follows directly from the two following lemmas.

**Lemma 7.19.** Let $(Z_{j,t})_{j \in \mathbb{N}^*, t \in \mathbb{Z}}$ be a family of $\mathcal{F}_{j,t}$-measurable random variables valued in $\mathbb{R}^{d'}$ for some positive integer $d'$ and satisfying (11). Suppose that either Assumption C1 or Assumption C'1 holds true.

Additionally, for a parameter $\theta$ belonging to a compact set $\Theta \subset \mathbb{R}^{d\theta}$, $n \ge 1$, let us consider a measurable function $m_{\theta}: Z \to \mathbb{R}$, and assume that

a) For all $z \in Z$, the application $\theta \mapsto m_{\theta}(z)$ is continuous.

b) For some $\delta > 0$, $\mathbb{E} \left( \sup_{\theta \in \Theta} |m_{\theta}(Z_{1,1})|^{1+\delta} \right) < +\infty$.

Then

$$ \sup_{\theta \in \Theta} \left| \frac{1}{nT} \sum_{t=1}^{T} \sum_{j=1}^{n} m_{\theta}(Z_{j,t}) - \mathbb{E}(m_{\theta}(Z_{1,1})) \right| \underset{n,T \to +\infty}{\longrightarrow} 0 \quad \text{a.s.} $$

*Proof of Lemma 7.19.* We set, for any $z \in Z$ and any $\Delta > 0$

$$ W_{\Delta}(z) := \sup \left\{ \left| m_{\theta+h}(z) - m_{\theta}(z) \right| : \theta \in \Theta, \theta+h \in \Theta, \|h\|_2 < \Delta \right\}. $$

One can easily show that

$$ \lim_{\Delta \to 0} \mathbb{E}(W_{\Delta}(Z_{1,1})) = 0. $$

Take any $\varepsilon > 0$, according to this latter equality, there exists $\delta > 0$ such that

$$ \Delta \le \delta \implies \mathbb{E}(W_{\Delta}(Z_{1,1})) < \varepsilon. $$
In addition, $\Theta \subset \cup_{\theta \in \Theta} B(\theta, \delta)$, where $B(\theta, \delta)$ denotes the open Euclidean ball of center $\theta$ and radius $\delta$. So, by compactness of $\Theta$, there exists a finite set $I$ such that $\theta \subset \cup_{i \in I} B(\theta_i, \delta)$.

For all $\theta \in \Theta$ and $n, T \in \mathbb{N}^*$, we denote

$$M_{n,T}(\theta) = \frac{1}{nT} \sum_{j=1}^{n} \sum_{t=1}^{T} m_{\theta}(Z_{j,t}) \quad \text{and} \quad M(\theta) = \mathbb{E}(m_{\theta}(Z_{1,1})).$$

For $\theta \in \Theta$, there exists $i \in I$ such that $\theta \in B(\theta_i, \delta)$ and for all $n, T \in \mathbb{N}^*$

$$ \begin{aligned} |M_{n,T}(\theta) - M(\theta)| &\le |M_{n,T}(\theta) - M_{n,T}(\theta_i)| + |M_{n,T}(\theta_i) - M(\theta_i)| + |M(\theta_i) - M(\theta)| \\ &\le \frac{1}{nT} \sum_{j=1}^{n} \sum_{t=1}^{T} |m_{\theta}(Z_{j,t}) - m_{\theta_i}(Z_{j,t})| + |M_{n,T}(\theta_i) - M(\theta_i)| + \mathbb{E}(W_{\delta}(Z_{1,1})) \\ &\le \frac{1}{nT} \sum_{j=1}^{n} \sum_{t=1}^{T} W_{\delta}(Z_{j,t}) + \max_{i \in I} |M_{n,T}(\theta_i) - M(\theta_i)| + \mathbb{E}(W_{\delta}(Z_{1,1})). \end{aligned} \quad (42) $$

▷ Using Assumption b) and either Proposition 4.1 or Proposition 4.2, we have for all $i \in I$

$$|M_{n,T}(\theta_i) - M(\theta_i)| \underset{n,T \to +\infty}{\longrightarrow} 0 \quad \text{a.s.} \quad (43)$$

We thus deduce that

$$\max_{i \in I} |M_{n,T}(\theta_i) - M(\theta_i)| \underset{n,T \to +\infty}{\longrightarrow} 0 \quad \text{a.s.}$$

▷ We then show that

$$\frac{1}{nT} \sum_{j=1}^{n} \sum_{t=1}^{T} W_{\delta}(Z_{j,t}) \underset{n,T \to +\infty}{\longrightarrow} \mathbb{E}(W_{\delta}(Z_{1,1})) \quad \text{a.s.} \quad (44)$$

We set here, for any $j \in \mathbb{N}^*$ and $t \in \mathbb{Z}$, $H_{j,t} = W_{\delta}(Z_{j,t})$. According to point b), for some $\kappa > 0$, we have

$$ \begin{aligned} \mathbb{E}(|H_{1,1}|^{1+\kappa}) &= \mathbb{E}(W_{\delta}(Z_{1,1})^{1+\kappa}) \\ &\le 2^{1+\kappa} \mathbb{E} \left( \sup_{\theta \in \Theta} |m_{\theta}(Z_{1,1})|^{1+\kappa} \right) \\ &< +\infty. \end{aligned} $$

Thus, Proposition 4.1 or Proposition 4.2 gives the convergence (44).

Inequality (42) then leads to

$$\sup_{\theta \in \Theta} |M_{n,T}(\theta) - M(\theta)| \le \frac{1}{nT} \sum_{j=1}^{n} \sum_{t=1}^{T} W_{\delta}(Z_{j,t}) + \max_{i \in I} |M_{n,T}(\theta_i) - M(\theta_i)| + \mathbb{E}(W_{\delta}(Z_{1,1}))$$

and by (43) and (44), we have almost surely

$$\sup_{\theta \in \Theta} |M_{n,T}(\theta) - M(\theta)| \le 2\mathbb{E}(W_{\delta}(Z_{1,1})) < 2\varepsilon$$

for $n$ and $T$ large enough, hence

$$\sup_{\theta \in \Theta} \left| \frac{1}{nT} \sum_{j=1}^{n} \sum_{t=1}^{T} m_{\theta}(Z_{j,t}) - \mathbb{E}(m_{\theta}(Z_{1,1})) \right| \underset{n,T \to +\infty}{\longrightarrow} 0 \quad \text{a.s.} \quad \square$$
**Lemma 7.20.** For some $i \in \{1,2\}$, let $(f_m)_{m \in \mathbb{N}^i}$ be a sequence of real-valued functions defined on a compact space $\Theta \subset \mathbb{R}^{d\Theta}$, such that $(f_m)_{m \in \mathbb{N}}$ converges uniformly toward a function $f$ on $\Theta$

$$ \sup_{\theta \in \Theta} |f_m(\theta) - f(\theta)| \underset{m \to +\infty}{\longrightarrow} 0. $$

Assume in addition that functions $f_m$, $m \in \mathbb{N}^i$ are continuous, and that $f$ admits a unique maximum in $\theta_0 \in \Theta$. Then

$$ \underset{\theta \in \Theta}{\operatorname{argmax}} f_m(\theta) \underset{m \to \infty}{\longrightarrow} \theta_0. $$

*Proof of Lemma 7.20.* Let $\|\cdot\|$ be an arbitrary norm on $\mathbb{R}^{d\Theta}$. For all $m \in \mathbb{N}^i$, let us denote $\hat{\theta}_m := \underset{\theta \in \Theta}{\operatorname{argmax}} f_m(\theta)$. Suppose that $\hat{\theta}_m$ does not converge to $\theta_0$. There then exists $\varepsilon > 0$ and a sequence $(m_k)_{k \in \mathbb{N}}$ valued in $\mathbb{N}^i$, with each coordinate going to $\infty$ and such that $\|\hat{\theta}_{m_k} - \theta_0\| > \varepsilon$. Since $\Theta$ is compact, up to an extraction, one can assume that for some $\theta^* \neq \theta_0 \in \Theta$,

$$ \hat{\theta}_{m_k} \underset{k \to \infty}{\longrightarrow} \theta^*. $$

We then have for all $k \in \mathbb{N}$

$$ |f_{m_k}(\hat{\theta}_{m_k}) - f(\theta^*)| \le |f_{m_k}(\hat{\theta}_{m_k}) - f(\hat{\theta}_{m_k})| + |f(\hat{\theta}_{m_k}) - f(\theta^*)|, $$

and the right side of this latter inequality converges toward 0 when $k \to \infty$ by uniform convergence and continuity. Thus

$$ f_{m_k}(\hat{\theta}_{m_k}) \underset{k \to \infty}{\longrightarrow} f(\theta^*). $$

Furthermore, for any $k \in \mathbb{N}$ we have

$$ f_{m_k}(\hat{\theta}_{m_k}) \ge f_{m_k}(\theta_0), $$

and considering the limits, it yields $f(\theta^*) \ge f(\theta_0)$. Since $\theta_0$ is the unique maximum of $f$ on $\Theta$, we deduce $\theta^* = \theta_0$. From this contradiction, we get the result. $\square$

## 7.10 Proof of Theorem 7.17

Once again, the proof of Theorem 7.17 is an adaptation of the proof of Theorem 7.10. We begin with a preliminary lemma.

**Lemma 7.21.** We have the following convergence in distribution

$$ \frac{1}{\sqrt{nT}} \sum_{t=1}^{T} \sum_{j=1}^{n} \nabla m_{\theta_0}(Z_{j,t}) \xrightarrow{\mathcal{L}} \mathcal{N}(0, \mathrm{Var}(\nabla m_{\theta_0}(Z_{1,1}))) $$

when $n, T \to +\infty$.

*Proof of Lemma 7.21.* Let us denote for all $n, T$ in $\mathbb{N}^*$

$$ S_{n,T} = \frac{1}{\sqrt{nT}} \sum_{t=1}^{T} \sum_{j=1}^{n} \nabla m_{\theta_0}(Z_{j,t}) \quad \text{and} \quad V = \mathrm{Var}(\nabla m_{\theta_0}(Z_{1,1})). $$

Note that $V$ is well defined according to assumption **P6**.

It is sufficient to show that for any increasing subsequences of integers $(n_k)_{k \ge 1}$ and $(T_k)_{k \ge 1}$, we have

$$ S_{n_k, T_k} \xrightarrow{\mathcal{L}} \mathcal{N}(0, V) \qquad (45) $$
when $k \to +\infty$.

Take indeed a function $g$ that is continuous and bounded, we want to show that

$$ \mathbb{E}(S_{n,T}) \underset{n,T \to +\infty}{\longrightarrow} \mathbb{E}(g(N)) $$

where $N \sim \mathcal{N}(0, V)$. If it was not the case, in the same spirit as in the proof of Theorem 7.15, we could build for some fixed $\varepsilon > 0$ two increasing subsequences of integers $(n_k)_{k \ge 1}$ and $(T_k)_{k \ge 1}$ such that

$$ \forall k \ge 1, |\mathbb{E}(g(S_{n_k, T_k})) - \mathbb{E}(g(N))| > \varepsilon, $$

which is in contradiction with (45).

We thus consider two increasing subsequences of integers $(n_k)_{k \ge 1}$ and $(T_k)_{k \ge 1}$ and denote for some $x \in \mathbb{R}^{d_\theta}$

$$ S_{n_k, T_k} = \sum_{t=1}^{T_k} \underbrace{\sum_{j=1}^{n_k} \frac{x' \nabla m_{\theta_0}(Z_{j,t})}{\sqrt{n_k T_k}}}_{V_{k,t}}, $$

so, using the Cramér-Wold device, we have to prove that

$$ \sum_{t=1}^{T_k} V_{k,t} \xrightarrow{\mathcal{L}} \mathcal{N}(0, x'Vx) $$

when $k \to +\infty$. To this end, we will apply Theorem 7.13 to the process $(W_{k,t})_{k,t \in \mathbb{N}^*}$ given for any $k, t \in \mathbb{N}^*$ by

$$ W_{k,t} = \begin{cases} V_{k,t} & \text{if } t \le T_k \\ 0 & \text{if } t > T_k \end{cases}. $$

The sum $\sum_{t=1}^{+\infty} W_{k,t}$ is necessarily finite for any $k \in \mathbb{N}^*$, and for each couple $(k, t)$, $W_{k,t}$ is in $L^2$, since it is the case for $V_{k,t}$ according to assumption **P10**. Thus we can define for any $k, t \in \mathbb{N}^*$

$$ \sigma_{k,t}^2 = \mathbb{E}(W_{k,t}^2 \mid \mathcal{H}_{t-1}) $$

where $\mathcal{H}_{t-1} = \bigvee_{j=1}^{+\infty} \mathcal{F}_{j,t-1}$. The family of $\sigma$-fields $(\mathcal{H}_t)_{t \in \mathbb{N}}$ is a filtration, and the sum $\sum_{t=1}^{+\infty} \sigma_{k,t}^2$ is also finite.

▷ From Assumption **P9**, the random sequence $(W_{k,t})_{t \in \mathbb{N}^*}$ is a martingale difference series adapted to the filtration $(\mathcal{H}_t)_{t \in \mathbb{N}^*}$.

▷ Let us then show that

$$ \sum_{t=1}^{+\infty} \sigma_{k,t}^2 \xrightarrow{\mathbb{P}} V $$

when $k \to \infty$. We actually have

$$ \begin{aligned} \sum_{t=1}^{\infty} \sigma_{k,t}^2 &= \sum_{t=1}^{T_k} \mathbb{E}(W_{k,t}^2 \mid \mathcal{H}_{t-1}) \\ &= \sum_{t=1}^{T_k} \sum_{j=1}^{n_k} \frac{1}{n_k T_k} \mathbb{E} \left( |x' \nabla m_{\theta_0}(Z_{j,t})|^2 \mid \mathcal{H}_{t-1} \right), \end{aligned} $$

since for $j \ne \ell$, $Z_{j,t}$ and $Z_{\ell,t}$ are conditionally independent given $\mathcal{H}_{t-1}$ (they depend on $\varepsilon_{j,t}$ and $\varepsilon_{\ell,t}$ and of some variables $\mathcal{H}_{t-1}$-measurable). Moreover setting $\kappa = |x' \nabla m_{\theta_0}|^2$, we have

$$ \tilde{Z}_{j,t} := \mathbb{E}(\kappa(Z_{j,t}) | \mathcal{H}_{t-1}) = \int \kappa \circ g (X_{j,t-1}, u, (\zeta_{t-s})_{s \ge 1}) \mathbb{P}_{\varepsilon_{1,1}}(du). $$