# Inferring the parameters of Taylor's law in ecology

Lionel Truquet (*), Joel E. Cohen (†), Paul Doukhan (‡)

August 30, 2024

## Abstract

Taylor's power law (TL) or fluctuation scaling has been verified empirically for the abundances of many species, human and non-human, and in many other fields including physics, meteorology, computer science, and finance. TL asserts that the variance is directly proportional to a power of the mean, exactly for population moments and, whether or not population moments exist, approximately for sample moments. In many papers, linear regression of log variance as a function of log mean is used to estimate TL's parameters. We provide some statistical guarantees with large-sample asymptotics for this kind of inference under general conditions, and we derive confidence intervals for the parameters. In many ecological applications, the means and variances are estimated over time or across space from arrays of abundance data collected at different locations and time points. When the ratio between the time-series length and the number of spatial points converges to a constant as both become large, the usual normalized statistics are asymptotically biased. We provide a bias correction to get correct confidence intervals. TL, widely studied in multiple sciences, is a source of challenging new statistical problems in a nonstationary spatiotemporal framework. We illustrate our results with both simulated and real data sets.

**Keywords:** asymptotic analysis of parameter estimates, biased confidence intervals, fluctuation scaling, space-time data, Taylor's law, variance function

## 1 Introduction

Taylor's law (TL) is an important family of relations between the mean and the variance of some probability distributions. TL is well approximated empirically in many data sets in the natural and social sciences, finance, and technology. Since Taylor (1961) publicized this relationship, which was discovered by others at least 20 years earlier, such as Bliss (1941), thousands of papers have studied TL for species abundance data in ecology and more generally to assess time-dependent or space-dependent changes of population distributions. Taylor (2019) and the references therein give a broad introduction to TL.

<footnote>*Lionel.TRUQUET@ensai.fr. Univ Rennes, Ensai, CNRS, CREST – UMR 9194, F-35000 Rennes, France</footnote>

<footnote>†cohen@rockefeller.edu. Laboratory of Populations, The Rockefeller University, New York, NY 10065, USA; also affiliated with Department of Statistics and Earth Institute, Columbia University, New York, NY 10027; and Department of Statistics, University of Chicago, Chicago, IL 60637, USA, ORCID 0000-0002-9746-6725</footnote>

<footnote>‡doukhan@cyu.fr. CY Cergy Paris University, AGM Mathematics, 2 Bd. Adolphe Chauvin, 95000 Cergy-Pontoise, AGM-UMR 8080</footnote>
The aim of this paper is to study the most basic version of this law, Taylor’s power law, as defined by (1) at the beginning of the next section. This law, which is based on a simple power relation between the variance and the mean of the marginal distribution, has an important interpretation in term of aggregation of "individuals" across time or space. This interpretation is probably one reason for its popularity in ecology. Cohen and Saitoh (2016) discuss the interpretation of the exponent (also called slope parameter) of Taylor’s power law. In scientific use and testing of TL, (1) is the most basic relation against many plausible alternative models. There are many extensions of TL such as the quadratic TL of Taylor et al. (1978) which can be more appropriate than (1), as shown, for instance, in Xu and Cohen (2021) for US county populations. Sometimes a cubic generalization of TL is required to describe data, as in Yang et al. (2022) for the study of mortality rates in Japan. And sometimes the variance function of the negative binomial distribution variance = a·mean+b·mean², as first proposed by Bartlett (1936), is superior to any form of TL, as shown in Cohen et al. (2017). Fitting TL (1) is a step that follows exploratory data analysis of multiple alternative models, each of which poses its own problems of estimating parameters.

Whatever the form of TL, estimation of unknown parameters roughly follows the same idea. We compute sample means and sample variances along one dimension (say, space) of a rectangular array of nonnegative data indexed by (say) space and time. This produces a pair (sample mean, sample variance) for each point in time. The sample mean and sample variance are assumed to be positive so their logarithms are well defined. The slope and intercept of log variance as a linear function of log mean are calculated, typically by ordinary least squares but sometimes by more robust methods (Yang et al., 2022) to optimize the fit of Taylor’s power law to these pairs of sample moments.

For example, in ecological species abundance data, the spatial TL, on which we will focus in the present paper, assumes that the time-varying mean and variance of the marginal distribution are linked at each time *t* by (1) and these population quantities are approximated by the corresponding sample moments computed from abundance data of a given species observed at different locations for each time *t*. After taking the logarithm of (1), this leads to the least-squares estimator given in Section 2.

A simple log-log plot (log variance on the vertical axis, log mean on the horizontal axis) is the classical graphical tool used by practitioners to assess whether TL describes the data. We show here that the confidence intervals of parameter estimates resulting from ordinary least squares may be seriously biased, and we show how to correct that bias. Section 6 illustrates with real data.

A wide statistical literature links power variance laws like TL with some classes of distributions (Tweedie, 1984; Bar-Lev and Enis, 1986; Bar-Lev and Stramer, 1987); those classes are compared in Bar-Lev (2020). Cohen and Huillet (2022) provide infinitely divisible examples of the Bar-Lev Enis type. However, formal statistical analyses of TL are rather scarce. Some methodological contributions include Clark and Perry (1994), James and Plank (2007), Jørgensen et al. (2011), Perry and Woiwod (1992) Sakoda et al. (2019), and Shi et al. (2017), among others. Beyond the previous papers which state that TL holds for interesting classes of distributions, these contributions mainly focus on inference procedures or computational issues and compare different approaches for estimating TL parameters in various situations such as sparsity, small samples, or data with covariates. None of these contributions discussed formally the large-sample behavior of least-squares or likelihood estimators, which is precisely the scope of the present paper.

We believe that the statistical analysis made here can guide the analysis of more complex versions of TL. Our initial intention was to study some detailed conditions that justify the commonly used method of log variance-log mean least-squares regression. However, our main assumptions and
results could be easily extended. Our approach to studying TL requires a careful analysis of a specific minimum of contrast estimator derived from some sample moments instead of the usual pair response/covariate used in standard regression analysis. Analyzing other methods such as likelihood estimation or more general models that use log-regression could be conducted using similar arguments.

This work can be useful for both quantitative ecologists and statisticians. Our analysis can help ecologists to measure uncertainty in the estimated parameters with precise confidence bands. It will prevent them from using standard but biased methods to compute confidence bands for TL parameters. By biased, we mean a non-centered asymptotic distribution for TL parameters. In Section 5, it is shown that in some cases, the Gaussian asymptotic distribution of least-squares estimates is not centered. For statisticians, our work reveals that a popular family of stochastic relations widely studied in several sciences is a source of interesting new theoretical problems, in particular for getting consistency properties or for constructing testing problems in a nonstationary spatiotemporal framework.

Many extensions of our results are possible, for example to the numerous variants of TL discussed above. In addition, this paper will focus only on distributions with light tails, i.e. when the probability distributions of interest have some finite moments of order $k > 2$. Though sufficient to prove rigorous consistency results for our estimates, such a restriction is not necessary for TL to occur. TL may hold for the sample mean and sample variance from stable laws with infinite variance or infinite mean, as shown in Brown et al. (2017), Brown et al. (2021), Cohen et al. (2020) or Cohen et al. (2022). However, even with our restriction to light-tailed distributions, ours is the first rigorous analysis of the asymptotic properties of the least-squares estimators of Taylor’s power law parameters.

The paper is organized as follows. Section 2 specifies our model and the least-squares estimators of TL parameters. Some numerical experiments in Section 3 assess the quality of these estimates, depending on various assumptions made about the process that generates the data. Section 4 gives consistency properties and Section 5 gives limiting distributions for our estimators. Section 6 illustrates the application of our results to real data. Section 7 summarizes our conclusions. All proofs of our results are postponed to Section 8. Section 9 contains some technical lemmas used in our proofs.

## 2 Least-squares estimators for Taylor’s law parameters

Let $(X_{j,t})_{1\le j\le n, 1\le t\le T}$ be an array of nonnegative random variables indexed by integers $j$, $t$. These variables can take integer values, for instance, when $X_{j,t}$ represents the number of individuals observed at time $t$ and at the location $j$, or more generally nonnegative real values, for instance, when $X_{j,t}$ represents the number of individuals observed at time $t$ divided by the area of site $j$, in which case $X_{j,t}$ represents a number of individuals per unit of area (population density).

We assume that the spatial variance $\sigma_t^2 := \text{Var}(X_{j,t})$ (that is, the variance of the $n$ elements in column $t$ of the array $(X_{j,t})$) and the spatial mean $\mu_t := \mathbb{E}(X_{j,t})$ (that is, the mean of the $n$ elements in column $t$ of the array $(X_{j,t})$) at time $t$ do not depend on the location $j$, are strictly positive, and satisfy

$$ \sigma_t^2 = \alpha \mu_t^\beta, \quad 1 \le t \le T, \qquad (1) $$

for a pair of parameters $(\alpha, \beta) \in (0, \infty) \times \mathbb{R}$. In most ecological applications, the estimate of $\beta$
is positive. A value $\beta > 1$ is sometimes interpreted as indicating aggregation or clumping of the organisms because a family of Poisson distributions in which the mean varies in time and remains constant over space obeys $\alpha = \beta = 1$. Values below 1 are uncommon but possible, for example, in models of a population in a Markovian environment, e.g., Cohen (2014a,b) and numerical simulations of some fisheries models gave $\beta < 0$ (Fujiwara and Cohen, 2015).

The population parameters ($\mu_t, \sigma_t^2$) are estimated by their empirical counterparts from the sample $X_{1,t}, \dots, X_{n,t}$. See below for details. This form of TL is known as a spatial TL because the empirical mean and the empirical variance depend on $n$ spatial locations at each time. We do not consider here other variants of TL, such as the temporal TL, in which the mean and the variance depend on $T$ time points at each spatial location, or the hierarchical spatial TL or the hierarchical temporal TL. For instance, the temporal TL assumes that the population moments $\sigma_j^2 = \text{Var}(X_{j,t})$ and $\mu_j = \mathbb{E}X_{j,t}$ depend only on $j$ and not on $t$. The corresponding power law is

$$ \sigma_j^2 = \alpha \mu_j^\beta, \quad 1 \le j \le n. $$

Asymptotic properties of the temporal TL can be obtained with similar arguments. This is why we mainly focus on the spatial law (1).

The main problem here is inference about the two parameters $\alpha$ and $\beta$ of (1). A natural way to proceed is to use a linear regression of the logarithm of the empirical spatial variance computed at each time $t$ as a function of the logarithm of the empirical spatial mean computed at each time $t$. Let

$$ \hat{\mu}_t = \frac{1}{n} \sum_{j=1}^{n} X_{j,t}, \quad \hat{\sigma}_t^2 = \frac{1}{n} \sum_{j=1}^{n} X_{j,t}^2 - \hat{\mu}_t^2. $$

While the variance can be also estimated from the unbiased estimator

$$ \hat{\sigma}_t^2 = \frac{1}{n-1} \sum_{j=1}^{n} (X_{j,t} - \hat{\mu}_t)^2, $$

our asymptotic results will not depend on the choice of a biased or unbiased version of the sample variance. Since (1) implies

$$ \log \sigma_t^2 = \log(\alpha) + \beta \log(\mu_t), \quad (2) $$

it is natural to estimate the pair $\theta = (\log(\alpha), \beta)'$ from the ordinary least-squares estimator

$$ \hat{\theta}_n = \arg\min_{\theta \in \mathbb{R}^2} \sum_{t=1}^T (\log \hat{\sigma}_t^2 - \theta_1 - \theta_2 \log \hat{\mu}_t)^2. $$

Explicitly,

$$ \hat{\theta}_n = \hat{D}_T^{-1} \hat{N}_T, \quad \text{where} \quad (3) $$

$$ \hat{D}_T = \left( \frac{1}{T} \sum_{t=1}^{T} \log \hat{\mu}_t \quad \frac{1}{T} \sum_{t=1}^{T} \log^2 \hat{\mu}_t \right), \quad \hat{N}_T = \left( \frac{1}{T} \sum_{t=1}^{T} \log \hat{\sigma}_t^2 \quad \log \hat{\mu}_t \cdot \log \hat{\sigma}_t^2 \right). \quad (4) $$

In the above formula, $\log \hat{\mu}_t$ or $\log \hat{\sigma}_t$ should be replaced by 0 if $\hat{\mu}_t = 0$ or $\hat{\sigma}_t^2 = 0$, respectively. For simplicity, we do not introduce an indicator function in the above formula. If we set

$$ D_T := \left( \frac{1}{T} \sum_{t=1}^{T} \log \mu_t \quad \frac{1}{T} \sum_{t=1}^{T} \log \mu_t \right), \quad N_T := \left( \frac{1}{T} \sum_{t=1}^{T} \log \sigma_t^2 \quad \log \mu_t \cdot \log \sigma_t^2 \right), $$
then analogously $\theta = D_T^{-1}N_T$, which makes it natural to introduce the above estimator. Throughout this paper, we will assume that the time series length $T$ depends on the number $n$ of spatial sites, i.e. $T := T_n$. We will study the asymptotic properties of $\hat{\theta} := \hat{\theta}_n$ in (3) as $n \to \infty$.

### 3 Numerical experiments

Before stating our main assumptions and results, we investigate the accuracy of the least-squares estimates discussed in the previous section through some numerical experiments. We will consider several scenarios for which the moment conditions on the marginal distribution are not the same. We also study a scenario with some zero-inflated data. Throughout this section, we consider a very simple time-varying parameter $\lambda_t = \exp(\cos(t))$, $1 \le t \le T$, to define our different scenarios. Moreover, for simplicity, we will assume that all the “observations” $\{X_{j,t}\}$ are independent across $j$ and $t$, i.e. there is no spatial or temporal dependence for the different random variables. To get a meaningful approximation of the population mean and population variance from empirical (or sample) counterparts, we assume that for any time point $t$, the random variables $X_{1,t}, \dots, X_{n,t}$ are identically distributed, i.e., all have the same distribution. We use DGP to abbreviate “data-generating process”. We use the relation $\sim$ always to mean that the random variables on both sides have the same probability distribution (¹).

1. **DGP1.** Poisson $X_{j,t} \sim P(\lambda_t)$. In this case, $\sigma_t^2 = \mu_t$ and $(\alpha, \beta) = (1, 1)$. All the exponential moments are finite here, i.e. $\mathbb{E}\exp(LX_{j,t}) < \infty$ for any $L > 0$.

2. **DGP2.** Squares of Gaussian $X_{j,t} \sim N(0, \lambda_t)^2$ (also known as $\chi^2$ with one degree of freedom). In this case, $\sigma_t^2 = 2\mu_t^2$ and $(\alpha, \beta) = (2, 2)$. Some exponential moments are infinite: $\mathbb{E}\exp(LX_{j,t}) < \infty$ if and only if $L\lambda_t < 1/2$.

3. **DGP3.** Mixture of Poisson $X_{j,t} \sim P(\lambda_t \times Z_{j,t} + \zeta_t)$, with $Z_{j,t} \sim \sqrt{3/4} \exp(N(0, 1)^2/8)$ and setting $v = \text{Var}(Z_{j,t})$,

$$ \zeta_t = \frac{1 + \sqrt{1 + 4v^2\lambda_t^2}}{2v} - \lambda_t. $$

In this case, $(\alpha, \beta) = (v, 2)$ and the fourth moment is infinite: $\mathbb{E}X_{j,t}^L < \infty$ if and only if $0 \le L < 4$.

<footnote n="1">¹Here we will never mean that deterministic functions on both sides are asymptotically equivalent as some quantity approaches a limit.</footnote>
<figure><img src="image_1.png" /><figcaption>(a) DGP1</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>(b) DGP2</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(c) DGP3</figcaption></figure>

<caption>Figure 1: Average RMSE values of the spatial TL for different values of *n* and *T*</caption>

Figure 1 represents the root mean squared error (RMSE) of $\hat{\beta}$ for the three scenarios and different values of the pair $(n, T)$. Large values of *n* and *T* seem to lead to a more accurate estimation of the slope parameter. However, the tail of the marginal distribution seems to affect this accuracy: the third scenario with infinite fourth moments leads to the largest estimation errors and large values of *T* can reduce the accuracy if the value of *n* is moderate. The difference in the accuracy of the
three DGPs can be seen from the differences in the vertical scales in the three panels of Figure 1. This observation was expected since, in theory, the moment conditions affect how closely empirical moments approximate population moments. It is not surprising to get smaller RMSE when the deviations $\hat{\mu}_t - \mu_t$ and $\hat{\sigma}_t^2 - \sigma_t^2$ can be made small.

We next evaluate how the number of 0 counts in the sample affects RMSE. Counts equal to 0 are very common in abundance data. For example, Figure 2 shows a real data set of fisheries abundances obtained from the ICES International Bottom Trawl Survey for the North Sea (NS-IBTS), which will be described in more detail in Section 5.

We now construct a new DGP with various proportions of 0 counts to mimic zero-inflated data. Let $W_{j,t}$, $1 \le j \le n$, $1 \le t \le T$, be Bernoulli random variables independently and identically distributed (iid) as $\mathcal{B}(p)$, $p \in (0, 1)$, where $p := \Pr\{W_{j,t} = 1\}$. Thus 0 values appear with probability $1-p$. Set $X_{j,t} = W_{j,t}X'_{j,t}$ where $X'_{j,t}$ are defined in **DGP2** and where $W_{j,t}$ are independent of $X'_{j,t}$. It may be shown that the spatial TL is still satisfied with $\beta = 2$. However, $\alpha = 3p - p^2$, which reduces to $\alpha = 2$ if $p = 1$.

Figure 3 shows that $p = 0.2$ gives less accurate results than $p = 0.5$. This observation is not surprising because smaller $p$ leads to more 0 values, and smaller values of the empirical means lead to larger absolute values for their logarithm and more variability is expected for the partial sums defining the log-regression estimates.

<figure><img src="image_4.png" /><figcaption>Figure 2: Time series for the species *Clupea harengus* on the subarea 31F2. CPUE is counts per unit effort.</figcaption></figure>
<figure><img src="image_5.png" /></figure>

$$ (a) X_{j,t} \sim \mathcal{N}(0, \lambda_t)^2 \times \mathcal{B}(0.2) $$

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

$$ (b) X_{j,t} \sim \mathcal{N}(0, \lambda_t)^2 \times \mathcal{B}(0.5) $$

<caption>Figure 3: Effect of zero-inflation on the RMSE of $\hat{\beta}$ for DGP2</caption>

# 4 Consistency properties

We now study the consistency of least-squares estimates. Consistency is defined to mean that $\lim_{n \to \infty} \hat{\theta} = \theta$ in probability (weak consistency) or almost surely (strong consistency).

We assume that the time horizon $T$ of the time series depends on $n$, i.e. $T = T_n$. If $T_+ := \sup_{n \ge 1} T_n < \infty$, consistency properties will roughly follow from the law of large numbers. It is more challenging to study the case $T_n \to \infty$, which is also useful in applications because the dimensions of the arrays of data are often balanced. Though time series with length converging to infinity will be the main focus, our results are also valid for finite $T_+$.

## 4.1 General results for nonnegative random variables

To derive the asymptotic properties of our estimates, we assume:

**A1** The time series $(X_{j,t})_{1 \le t \le T}$, $1 \le j \le n$, are independent and identically distributed.

**A2** There exists $k > 2$ such that $\sup_{t \ge 1} \mathbb{E}|X_{1,t}|^k < \infty$.
**A3** There exists some $\delta > 0$ such that

$$p(\delta) := \inf_{t \ge 1} \mathbb{P}(X_{1,t} > \delta) > 0 \quad \text{and} \quad q(\delta) = \inf_{t \ge 1} \mathbb{P}(|X_{1,t} - X_{2,t}| > \delta) > 0.$$

**A4** There exists $c > 0$ such that for any positive integer $T$,

$$\frac{1}{T} \sum_{t=1}^{T} (\log \mu_t)^2 - \left( \frac{1}{T} \sum_{t=1}^{T} \log \mu_t \right)^2 \ge c.$$

Let us comment on our assumptions.

1. Assumption **A1** does not require any specific dependence property across time. Independence across sites $j$ is needed only for technical reasons. However, spatial correlations can be observed in data, as for instance in a study of Hokkaido voles (Cohen and Saitoh, 2016). In theory, it may be possible to remove the assumption of independence across sites $j$ by assuming a stationarity condition and some weak dependence property between the different sites $j$. However, getting consistency properties under such conditions would require additional assumptions, in particular to get the law of large numbers when the locations $j$ are irregularly spaced. We prefer to keep the setup simple.

2. The technical moment assumption **A2** is needed to prove limit results as in Theorem 1.

3. Assumption **A3** roughly means that Dirac masses cannot be cluster points of the marginal distribution, since **A3** is not satisfied if and only if at least one of the two sequences $(X_{1,t})_{t \ge 1}$ or $(X_{1,t} - X_{2,t})_{t \ge 1}$ has a subsequence converging to 0 in probability. Excluding Dirac masses is important to avoid small values for the population mean and the population variance as well as for their empirical versions and then large negative values for their logarithm.

4. Assumption **A4** ensures that the vectors $(\log \mu_t)_{1 \le t \le T}$ vary sufficiently in time. Such a condition is needed to avoid some pathological behaviors such as the convergence of $\mu_t$ when $t \to \infty$. Without this condition, the matrix $\hat{D}_T$ can be asymptotically singular. For a standard linear regression model $y_t = x'_t\beta + \varepsilon_t$, our condition is comparable to the existence of a lower bound for the determinant of the matrices $T^{-1} \sum_{t=1}^{T} x_t x'_t$ classical in statistical learning.

**Theorem 1.** Given Assumptions **A1-A2-A3-A4**, there exists $\gamma > 0$ such that if $T_n = O(\exp(\gamma n))$, then $\lim_{n \to \infty} \hat{\theta} = \theta$ in probability. If Assumption **A2** holds with $k > 4$, then $\lim_{n \to \infty} \hat{\theta} = \theta$ almost surely (a.s.).

Theorem 1 requires a limited growth of the time series length $T$ with respect to the number $n$ of sites. Another assumption concerning the marginal distribution of $X_{1,t}$ avoids such a restriction and leads to consistency. Define $\log_-(x) := -\log \min(x, 1)$ for any $x > 0$ and $\log_-(x) := (\log_-(x))^q$ for any $q > 0$.

**A5** There exists $q > 2$ such that

$$\sup_{t \ge 1} \mathbb{E}1_{X_{1,t} > 0} \cdot \log_-' X_{1,t} < \infty, \quad \sup_{t \ge 1} \mathbb{E}1_{X_{1,t} \ne X_{2,t}} \cdot \log_-' |X_{1,t} - X_{2,t}| < \infty.$$

**Theorem 2.** Given Assumptions **A1-A2-A3-A4-A5**, $\hat{\theta}_n$ is weakly consistent if $k > 2$ and strongly consistent if $k > 4$ in **A2**.
## 4.2 A more specific result for integer-valued random variables

For integer-valued random variables, the previous assumptions can be simplified. In particular, it is straightforward to show that Assumption A5 is automatically satisfied. Assumption A3 becomes:

**A3'** We have $\inf_{t \ge 1} \mathbb{P}(X_{1,t} > 0) > 0$ and $\inf_{t \ge 1} \mathbb{P}(X_{1,t} \neq X_{2,t}) > 0$.

Then:

**Corollary 1.** Given Assumptions A1-A2-A3'-A4, $\hat{\theta}_n$ is weakly consistent if $k > 2$ and strongly consistent if $k > 4$ in A2.

# 5 Asymptotic normality of parameter estimates

We now derive asymptotic normality of our log-regression OLS estimators. To keep the setup as simple as possible, we will also assume independence in the time direction, i.e. for any $j \ge 1$, the $X_{j,1}, X_{j,2}, \dots$ are assumed to be mutually independent. This assumption can be relaxed but this simple setup suffices to give insight into the weak convergence of such estimates.

Define $Y_{j,t} := (X_{j,t}, X_{j,t}^2)'$ and $m_t := \mathbb{E}Y_{1,t}$ for $1 \le j \le n$ and $1 \le t \le T$. Define the empirical version of $m_t$ by

$$ \hat{m}_t := \frac{1}{n} \sum_{j=1}^{n} Y_{j,t}. $$

The two components of $m_t$ and $\hat{m}_t$ will be denoted respectively by $m_{1,t}, m_{2,t}$ and $\hat{m}_{1,t}, \hat{m}_{2,t}$. In terms of our previous notations, we have $m_{1,t} = \mu_t$ and $\hat{m}_{1,t} = \hat{\mu}_t$. In what follows, we will use both notations, depending on the context.

We also use a standard notation for weak convergence. When $(Z_n)_{n \ge 1}$ is a sequence of random vectors of $\mathbb{R}^2$ and $\nu$ is a probability measure on the Borel $\sigma$-field of $\mathbb{R}^2$, notation $Z_n \Rightarrow \nu$ means that $Z_n$ converges weakly to $\nu$.

In the following analysis, we are thinking of a function $\phi_\theta$ (with parameter $\theta \in \mathbb{R}^2$) from domain $(0, \infty)^2$ into range $\mathbb{R}^2$. In the 2-dimensional argument $(x, y)$ of $\phi_\theta$, $x$ represents the empirical or expected spatial mean of the observations $X_{j,t}$ and $y$ represents the empirical or expected spatial mean of the squared observations $X_{j,t}^2$. With either empirical or expected mean, $y - x^2$ represents the variance, i.e., the mean squared observation minus the squared mean observation. The two parameters $\theta_1 = \log(\alpha)$, $\theta_2 = \beta$ of $\phi_\theta$ correspond to the intercept and the slope of the log-log form (2) of TL. The first component $\phi_{1,\theta}(x, y) := \log(y - x^2) - \theta_1 - \theta_2 \log(x)$ of $\phi_\theta$ is the deviation from (2), conditional on the two parameter values. It corresponds to the first component of the normal equations of ordinary least squares estimation, (3) and (4). The second component $\phi_{2,\theta}(x, y) := \log(x) \log(y - x^2) - \theta_1 \log(x) - \theta_2(\log(x))^2$ similarly corresponds to the second component of the normal equations of ordinary least squares.

Setting $\theta = (\log(\alpha), \beta)'$, the estimation error for the log-regression parameters has the following decomposition:

$$ \begin{aligned} \hat{\theta} - \theta &= \hat{D}_T^{-1} \hat{N}_T - D_T^{-1} N_T \\ &= \hat{D}_T^{-1} \left( \hat{N}_T - N_T + (D_T - \hat{D}_T) \theta \right). \end{aligned} $$
This decomposition implies that

$$ U_T := \hat{D}_T(\hat{\theta} - \theta) = \frac{1}{T} \sum_{t=1}^{T} (\phi_{\theta}(\hat{m}_t) - \phi_{\theta}(m_t)), \qquad (5) $$

where

$$ \phi_{\theta}(x, y) := \begin{pmatrix} \phi_{1,\theta}(x, y) \\ \phi_{2,\theta}(x, y) \end{pmatrix}, \quad x > 0, y > 0. $$

The idea is simply to make a Taylor expansion of the function $\phi_{\theta}$ with respect to its two arguments representing the spatial means of the two samples $(X_{j,t})_{1 \le j \le n}$ and $(X_{j,t}^2)_{1 \le j \le n}$ respectively, computed at any time point $1 \le t \le T$. We will then approximate $U_T$ with its first-order Taylor expansion

$$ V_T = \frac{1}{T} \sum_{t=1}^{T} J_{\phi_{\theta}}(m_t) \cdot (\hat{m}_t - m_t), \qquad (6) $$

where $J_{\phi_{\theta}}(x, y)$ denotes the Jacobian matrix of $\phi_{\theta}$ at point $(x, y)$. Intuitively, the convergence rate for $V_T$ is $\sqrt{nT}$ since $V_T$ contains a double summation with $n \times T$ centered terms. We then define

$$ \Gamma_T := \mathrm{Var}(\sqrt{nT}V_T) = \frac{1}{T} \sum_{t=1}^{T} J_{\phi_{\theta}}(m_t) \cdot \mathrm{Var}(Y_{1,t}) \cdot J_{\phi_{\theta}}(m_t)' \qquad (7) $$

and $C_T := \Gamma_T^{-1/2}$. It is expected that $\sqrt{nT}C_TV_T \Rightarrow \mathcal{N}_2(0, I_2)$. The remaining step is to investigate the asymptotic behavior of $\sqrt{nT}(U_T - V_T)$. When $T/n \to 0$, one can show that under the assumptions of Theorem 1, $\sqrt{nT}(U_T - V_T) = o_P(1)$. However, when $T/n \to s > 0$, this is not the case, but under the assumptions of Theorem 1, one can show that

$$ \sqrt{nT}(U_T - V_T) = \sqrt{\frac{n}{T}} \sum_{t=1}^{T} \left( \begin{pmatrix} (\hat{m}_t - m_t)' \nabla^{(2)} \phi_{1,\theta}(m_t) (\hat{m}_t - m_t) \\ (\hat{m}_t - m_t)' \nabla^{(2)} \phi_{2,\theta}(m_t) (\hat{m}_t - m_t) \end{pmatrix} + o_P(1), \right) $$

where for a twice-differentiable function $f: \mathbb{R}^2 \to \mathbb{R}$, we denote by $\nabla^{(2)} f(m)$ the Hessian matrix of $f$ at point $m$. The first right-hand term in the previous equality does not vanish asymptotically and has an expectation of order $\sqrt{T/n}$. Hence this non-negligible bias has to be included for deriving the asymptotic distribution of $\hat{\theta}$.

**Theorem 3.** Let $\lambda_-(\Gamma_T)$ denote the smallest eigenvalue of $\Gamma_T$ defined in (7). Suppose that Assumptions A1-A2-A3-A4 hold with $k > 4$ in A2 and for any integer $j \ge 1$, $(X_{j,t})_{t \ge 1}$ is a sequence of independent random variables. Assume furthermore that $\lim_{n} \lambda_-(\Gamma_T) > 0$.

1. If $T_n/n \to 0$, then

$$ \sqrt{nT}C_T \hat{D}_T (\hat{\theta} - \theta) \Rightarrow \mathcal{N}_2(0, I_2). $$

2. If $T_n/n = O(1)$, then

$$ \sqrt{nT}C_T \hat{D}_T (\hat{\theta} - \theta) - \frac{1}{2} \sqrt{\frac{T}{n}} C_T E_T \Rightarrow \mathcal{N}_2(0, I_2), $$
where $E_T := \frac{1}{T} \sum_{t=1}^{T} G_t$ and

$$G_t := \left( \mathrm{Tr}(\nabla^{(2)} \phi_{1,\theta}(m_t) \mathrm{Var}(Y_{1,t})), \mathrm{Tr}(\nabla^{(2)} \phi_{2,\theta}(m_t) \mathrm{Var}(Y_{1,t})) \right)'$$

We comment on these results.

1. The existence of a lower bound for the eigenvalues of $\Gamma_T$ is always difficult to check for a specific sequence of marginal distributions, though quite natural. We provide sufficient conditions as well as a way to check it for Poisson random variables in Lemma 6 in the Appendix.

2. When $T_n/n \to k > 0$ as $n \to \infty$, the normalized statistics $\sqrt{nT}C_T\hat{D}_T(\hat{\theta} - \theta)$ in Theorem 3.1 become asymptotically biased. Neither the matrix $C_T$ nor the vector $E_T$ converges in general when $n \to \infty$ but these sequences remain bounded under our assumptions. This non-centering of the asymptotic distribution is atypical for regular statistical models and is due here to the fact that when the array of data remains balanced, there are not enough replicates of the time series to make the second-order term in the Taylor expansion of the estimates negligible. The asymptotics when $T/n = O(1)$ (Theorem 3.2) also hold when $T/n = o(1)$, since then $\sqrt{T/n}C_T E_T = o(1)$ becomes negligible.

3. Figure 4 shows the effect of neglecting the bias for simulated Poisson data when the number of sites equals the number of time points. Without bias correction, the probability distribution of the two coordinates of the normalized statistics $\sqrt{nT}\hat{C}_T\hat{D}_T(\hat{\theta} - \theta)$ are not well approximated by that of a standard Gaussian.

4. To use these results in practice for constructing tests or confidence bands, it is necessary to estimate $\Gamma_T$ and $E_T$ from the data. This can be done with no additional assumptions:

**Corollary 2.** Suppose that all the assumptions of Theorem 3 hold true with either the asymptotics $T/n = o(1)$ or $T/n = O(1)$. Set $\Sigma_t := \mathrm{Var}(Y_{1,t})$ and

$$\hat{\Sigma}_t := \frac{1}{n} \sum_{j=1}^{n} Y_{j,t} Y'_{j,t} - \bar{Y}_{n,t} \bar{Y}'_{n,t}, \quad \bar{Y}_{n,t} := \frac{1}{n} \sum_{j=1}^{n} Y_{j,t},$$

$$\hat{\Gamma}_T := \frac{1}{T} \sum_{t=1}^{T} J_{\phi_{\hat{\theta}}}(\hat{m}_t) \hat{\Sigma}_t J_{\phi_{\hat{\theta}}}(\hat{m}_t)', \quad \hat{C}_T := \hat{\Gamma}_T^{-1/2},$$

$$\hat{E}_T := \frac{1}{T} \sum_{t=1}^{T} \hat{G}_t,$$

$$\hat{G}_t := \left( \mathrm{Tr}(\nabla^{(2)} \phi_{1,\hat{\theta}}(\hat{m}_t) \hat{\Sigma}_t), \mathrm{Tr}(\nabla^{(2)} \phi_{2,\hat{\theta}}(\hat{m}_t) \hat{\Sigma}_t) \right)'$$

Then

$$\sqrt{nT}\hat{C}_T\hat{D}_T(\hat{\theta} - \theta) - \frac{1}{2}\sqrt{\frac{T}{n}}\hat{C}_T\hat{E}_T \Rightarrow \mathcal{N}_2(0, I_2). \qquad (8)$$
Define $H_n := \hat{D}_T^{-1}\hat{\Gamma}_T\hat{D}_T^{-1}$ and

$$M_n := \begin{cases} 0 & \text{if } T/n = o(1), \\ \sqrt{\frac{T}{n}}\hat{D}_T^{-1}\hat{E}_T & \text{if } T/n = O(1). \end{cases}$$

Let $q_β$ denote the quantile of order $β$ of the standard Gaussian distribution. Then the interval

$$I_{i,n} = \left[ \hat{\theta}_i - \frac{\sqrt{H_n(i,i)}q_{1-\alpha/2} + M_n(i)}{\sqrt{nT}}, \hat{\theta}_i + \frac{\sqrt{H_n(i,i)}q_{1-\alpha/2} - M_n(i)}{\sqrt{nT}} \right]$$

is a confidence interval with asymptotic level $1 - \alpha$ for $i = 1, 2$.

In practice, it is difficult to know which regime is the most appropriate. However, if $T/n = o(1)$, we have $M_n = o(1)$, under our assumptions. Then the bias correction can still be used in this case, though it is asymptotically negligible. Moreover, even when $T/n$ seems to be small, the bias correction $M_n$ depends on some empirical averages and can be non-negligible. This is why we recommend to use such a bias correction whatever the context, or at least to compute both confidence intervals and to see how they differ from each other.
<figure><img src="image_7.png" /><figcaption>Figure 4: For simulated Poisson variates $X_{j,t} \sim P(\lambda_t)$ with $n = T = 100$ and $\lambda_t = 3 + \cos(t)$, quantile-quantile plots (top) for the two components of statistics $\sqrt{nT}\hat{C}_T\hat{D}_T(\hat{\theta} - \theta)$ and (bottom) for the two components of statistics $\sqrt{nT}\hat{C}_T\hat{D}_T(\hat{\theta} - \theta) - \frac{1}{2}\sqrt{T/n}\hat{C}_T\hat{E}_T$.</figcaption></figure>

# 6 Example: International Bottom Trawl Survey for the North Sea

We analyze some data from the International Council for the Exploration of the Sea (ICES) International Bottom Trawl Survey for the North Sea (NS-IBTS). The survey estimates the abundance, measured as Catch Per Unit Effort (CPUE), of each fish and some invertebrate species binned into length classes within spatial rectangles of 0.5° latitude by 1.0° longitude, with an approximate area of 30 km², known as subareas. Data included in this study are limited to the years 1977 through 2015. The dataset used here contains only the observations of the nine standard species: *Clupea harengus* (species 1), *Gadus morhua* (species 2), *Melanogrammus aeglefinus* (species 3),
*Merlangius merlangus* (species 4), *Pleuronectes platessa* (species 5), *Pollachius virens* (species 6), *Scomber scombrus* (species 7), *Sprattus sprattus* (species 8) and *Trisopterus esmarkii* (species 9). The CPUE is given by the number of individuals caught per cruise in 30 minutes of fishing, where only individuals in certain categories of length (measured in mm) were included. Further, data were limited to quarter one only, when survey cruises typically occur from mid January to mid February. It then excludes seasonal differences in abundance. For each of the nine species, the abundance is recorded at *n* = 195 locations and across *T* = 39 years, with one observation per year. See Cobain *et al.* (2019) and the references therein for more details about this survey.

Figure 5 gives log-log plots (using the natural logarithm) of the empirical mean and the empirical variance of the CPUE of the nine species.

Table 1 gives the OLS estimate of the slope parameter as well as the 95% confidence bands obtained from standard linear regression methods. In contrast, Table 2 provides confidence bands computed from our asymptotic results, with and without bias correction, for the same confidence level. We do not consider species 7 here because more than 75% of these data were zeros and we had numerical difficulties evaluating the quantities in Corollary 2. Also, we divided all the abundances by a factor (the mean value of the abundance of the species, computed over the $n \times T$ observations) to avoid numerical issues due to many large values in the time series. This operation does not change the value of the estimate of the slope parameter $\beta$ because TL is still valid with the same $\beta$, although it changes the value of the intercept parameter $\alpha$.

The confidence intervals given in Table 1 are wider. This is not surprising since their size is proportional to $1/\sqrt{T}$ while the correct size is proportional to $1/\sqrt{nT}$. As a consequence, after this correction of the width of the confidence interval, some values of $\beta$ are no longer admissible at the 95% level, for instance, the value $\beta = 2$ for Species 4. Recall that the two intervals with or without bias correction have exactly the same length. The confidence interval obtained with bias correction is a simple shift, to the left or to the right, of that obtained without bias correction and it is then no longer centered around the estimate $\hat{\beta}$.

With or without bias correction, the confidence intervals of four species, 1, 5, 8, 9, include the value $\beta = 2$; three intervals lie entirely above $\beta = 2$; and two intervals lie entirely below $\beta = 2$. Thus the hypothesis $\beta = 2$ is not rejected for only four species, under the important, probably influential assumption, for each species, that its time series are independently and identically distributed at different spatial locations. Our methods (and most other applications of least-squares regression to estimate the parameters of TL) ignore the possibility that exogenous environmental fluctuations like the North Atlantic Oscillation influence all locations simultaneously, as well as the possibility that spatially adjacent or spatially remote populations of a species interact.
<figure><img src="image_8.png" /><figcaption>(a) Species 1</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) Species 2</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(c) Species 3</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(d) Species 4</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(e) Species 5</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(f) Species 6</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(g) Species 7</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(h) Species 8</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(i) Species 9</figcaption></figure>

Figure 5: Plots of the natural logarithm of sample variances versus the natural logarithm of sample means of the CPUE of the nine species. The ordinary least-squares regression lines are red.
<table><caption>Table 1: Point estimate by ordinary least squares of the slope parameter $\beta$ for the nine species and lower and upper bound of the 95% confidence interval obtained with the conventional linear regression</caption><thead><tr><th>Species</th><th>OLS β</th><th>lower bound</th><th>upper bound</th></tr></thead><tbody><tr><td>1</td><td>2.009793</td><td>1.780779</td><td>2.238807</td></tr><tr><td>2</td><td>2.390582</td><td>2.1058999</td><td>2.675264</td></tr><tr><td>3</td><td>1.530209</td><td>1.221486</td><td>1.838933</td></tr><tr><td>4</td><td>1.696704</td><td>1.301093</td><td>2.092315</td></tr><tr><td>5</td><td>2.106150</td><td>1.780474</td><td>2.431825</td></tr><tr><td>6</td><td>2.205332</td><td>2.041850</td><td>2.368813</td></tr><tr><td>7</td><td>2.09691</td><td>2.015611</td><td>2.178211</td></tr><tr><td>8</td><td>2.084444</td><td>1.779885</td><td>2.389003</td></tr><tr><td>9</td><td>1.893228</td><td>1.643332</td><td>2.143124</td></tr></tbody></table>

<table><caption>Table 2: Point estimate by ordinary least squares of the slope parameter $\beta$ for eight of the nine species and lower and upper bound of the 95% confidence interval ($n = 195$, $T = 39$) computed with bias correction (C) and without bias correction (NC)</caption><thead><tr><th>Species</th><th>OLS β</th><th>lower bound NC</th><th>upper bound NC</th><th>lower bound C</th><th>upper bound C</th></tr></thead><tbody><tr><td>1</td><td>2.009793</td><td>1.886288</td><td>2.133297</td><td>1.918103</td><td>2.165112</td></tr><tr><td>2</td><td>2.390582</td><td>2.228910</td><td>2.552254</td><td>2.286405</td><td>2.609748</td></tr><tr><td>3</td><td>1.530209</td><td>1.379527</td><td>1.680892</td><td>1.346237</td><td>1.647602</td></tr><tr><td>4</td><td>1.696704</td><td>1.429183</td><td>1.964225</td><td>1.356336</td><td>1.891378</td></tr><tr><td>5</td><td>2.106150</td><td>1.926557</td><td>2.285742</td><td>1.965488</td><td>2.324674</td></tr><tr><td>6</td><td>2.205332</td><td>2.073231</td><td>2.337432</td><td>2.191456</td><td>2.455656</td></tr><tr><td>8</td><td>2.084444</td><td>1.864855</td><td>2.304033</td><td>1.983217</td><td>2.422395</td></tr><tr><td>9</td><td>1.893228</td><td>1.735085</td><td>2.051371</td><td>1.770439</td><td>2.086725</td></tr></tbody></table>

Visual inspection of Figure 5 gives no suggestion, for any species, that log variance is systematically non-linearly related to log mean. In an exploratory analysis, we fitted a quadratic spatial TL, namely, log spatial variance = $\log(a) + b \times (\log spatial mean) + c \times (\log spatial mean)^2$, by ordinary least squares as implemented in the software R. The null hypothesis that $c = 0$ was not rejected for any species, but the confidence intervals for $c$ given by the software are not credible in the context of this application. Further theoretical developments are required to assess rigorously whether the data reject the null hypothesis of a linear relation between log spatial mean and log spatial variance. For now, we accept the visual evidence that TL, linear on log-log coordinates, is a plausible model of the variance function of these data.

The present paper does not analyze stochastic dependence between observations. Considering dependence will lead to a more delicate estimation of the normalizing matrix $C_T$ in Theorem 3 and it would probably require large data sets to get accurate results.

To decide if this data set is compatible with the assumption of temporal independence under a nonstationary framework, we compute the residuals $\hat{Z}_{j,t} = \frac{X_{j,t} - \hat{\mu}_t}{\hat{\sigma}_t}$ which behave more like a
stationary time series since their two first moments approximate 0 and 1. The arctanh function has been applied to all the correlation coefficients as a variance-stabilizing function. For each species, we then compute the three first autocorrelation coefficients $\rho(1)$, $\rho(2)$, $\rho(3)$ with lags 1, 2, 3 of each time series $(\hat{Z}_{j,t})_{1\le t\le T}$ for $j = 1, \dots, n$. Van der Vaart (2000, Section 3.2) defines the corresponding confidence interval of the autocorrelation coefficients, based on this variance-stabilizing function.

Table 3 gives the percentage of confidence bands that include 0 for each species. The confidence bands were computed under the assumption of independence across $t$ of $Z_{j,t} = \frac{X_{j,t}-\mu_t}{\sigma_t}$ and by neglecting the sampling error in estimating the local mean and local variance. Because $n = 195$ is much larger than $T = 39$ for this example, one might hope this sampling error would be negligible. However, the results have to be taken with caution since $T$ is not very large. For species 6 to 9, the assumption of no correlation across time seems to be perfectly plausible. Species 3 seems to be the most problematic.

In Table 4, we also apply the same procedure to test spatial correlation. More precisely, for each of the nine species, we compute the sample or empirical correlation $r$ between $Z_{j,t}$ and $Z_{\ell,t}$ for two different locations $j$ and $\ell$ using the two samples $(Z_{j,t})_{1\le t\le T}$ and $(Z_{\ell,t})_{1\le t\le T}$. We report the percentage of pairs $(j, \ell)$ for which $r = 0$ is not rejected. As for temporal dependence, we note that for some species, such as species 6, it is difficult to exclude a spatial correlation whereas for some other, spatial correlation seems not to be very important.

Adding an assumption of temporal dependence and or spatial dependence could be an interesting extension of this work. On the other hand, our empirical results suggest that these independence assumptions are a reasonable first approximation for this data set.

We have given a first rigorous basis for estimating the TL parameters and we have shown that our methods are applicable and illuminating for a real data set. When extending our work to allow spatial and temporal dependence, a similar theoretical analysis can be conducted using our proof procedure. The asymptotics of the ratio $T/n$ will be also of crucial importance in the statement of the corresponding results.

<table><caption>Table 3: Percentage of confidence intervals containing the point $\rho(i) = 0$ for serial autocorrelation with lag $i = 1, 2, 3$. The confidence intervals were computed for the $n = 195$ time series using the arctanh variance-stabilizing function.</caption><thead><tr><th>Species</th><th>ρ(1)</th><th>ρ(2)</th><th>ρ(3)</th></tr></thead><tbody><tr><td>1</td><td>77.4</td><td>80.0</td><td>99.0</td></tr><tr><td>2</td><td>75.4</td><td>89.2</td><td>91.3</td></tr><tr><td>3</td><td>48.2</td><td>51.3</td><td>52.8</td></tr><tr><td>4</td><td>70.0</td><td>79.0</td><td>82.0</td></tr><tr><td>5</td><td>81.5</td><td>85.6</td><td>90.3</td></tr><tr><td>6</td><td>97.9</td><td>98.0</td><td>100.0</td></tr><tr><td>7</td><td>95.9</td><td>98.5</td><td>99.0</td></tr><tr><td>8</td><td>97.4</td><td>97.4</td><td>96.9</td></tr><tr><td>9</td><td>95.4</td><td>95.4</td><td>96.9</td></tr></tbody></table>
<table><caption>Table 4: Percentage of confidence intervals containing the point $r = 0$ for correlation of the abundance between different locations. The confidence intervals were computed using the arctanh variance-stabilizing function.</caption><thead><tr><th>Species</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th></tr></thead><tbody><tr><th>Percentage</th><td>83</td><td>84.5</td><td>60.4</td><td>83.8</td><td>63.4</td><td>49.4</td><td>68.1</td><td>66.7</td><td>66</td></tr></tbody></table>

# 7 Conclusions and perspectives

In this work, we introduced precise assumptions on the data-generating process to study the asymptotic properties of least-squares estimators of Taylor's power law parameters. Theorems 1, 2, and 3 are our main new results. Our assumptions on the marginal distributions are general enough to cover many situations encountered in practice, excluding data with fat tails and spatial and/or temporal dependence.

To ecologists interested in using our results to control the uncertainty in estimates of TL parameters, we offer the following guidance.

* Always apply the confidence intervals corresponding to the bias correction given in Corollary 2, with $M_n$ defined as if $T/n = O(1)$. This correction is crucial when $T$ and $n$ are of the same magnitude. Moreover, even if the ratio $T/n$ seems to be quite small, it is difficult to discuss a particular threshold to know if the asymptotic regime $T/n = o(1)$ is appropriate. The formula with bias correction always works asymptotically when $T/n = o(1)$. Even for finite sample sizes, the confidence bands should not be too much affected.

* The most important distortion of the confidence intervals could come from spatial or temporal dependence that we ignored in this work. We recommend checking whether the independence assumption is compatible with the data set, as we did, for instance, in our real data application. If an important deviation from independence is observed, further research is needed. Our specific formulas for confidence intervals may or may not prove to be applicable to data with demonstrable spatial and/or temporal dependence, but we believe our approach to modeling and analysis of TL parameter estimation generalizes even to data with dependence.

Some hints for various extensions are given below. When data have fat tails, the limiting distribution of least-squares estimators will probably no longer be Gaussian. Even with thin tails, our contribution appears to be the first to derive an asymptotic distribution for such estimates, which appear in numerous ecological papers without theoretical analysis.

Various extensions could also be investigated. For instance, one can study other models

$$ \log \sigma_t^2 = f_{\theta}(\log \mu_t) \qquad (9) $$

where $f_{\theta}$ is an unknown function depending on a finite number of parameters, such as the quadratic and the cubic TL mentioned in the introduction. Another extension could use a statistical learning point of view by studying more general least-squares problems of the form

$$ \hat{f} = \arg \min_{f \in \mathcal{F}} \sum_{t=1}^{T} (\log \hat{\sigma}_t^2 - f(\log \hat{\mu}_t))^2, \qquad (10) $$
where $\mathcal{F}$ is an infinite-dimensional class of functions. Of course, we have to take care that whatever the regularity (e.g. differentiability) imposed on $\mathcal{F}$, there will be an infinite number of functions mapping the sequence $(\log \mu_t)_t$ onto the sequence $(\log \sigma_t^2)_t$. This non-parametric extension (10) will differ from the usual standard non-parametric estimation problems in statistics. Testing a specific parametric model of the form (9) could be also interesting.

Another open question related to Assumption **A1** is to check whether time and space $(t$ and $j$) could be exchanged. Non-stationarity in time seems mandatory to get identifiable parameters $\alpha, \beta$ in (1). Estimating the normalizing quantities could be more difficult in this nonstationary context. Taking account of the stochastic dependence, spatial or temporal or both, between the observations is also a technical challenge. Many ways to extend Theorem 3 in the context of dependence have been published; for instance, Dedecker et al. (2007) provide many helpful asymptotic results.

Extensions to dependent random variables are plausible in different settings. For example, de la Pena et al. (2022) suggest a variant definition of TL together with a simple extension of the TL under mixing type conditions. Cohen et al. (2020) consider heavy-tailed processes. Doukhan et al. (2002) give a systematic approach to a weaker condition such as long-range dependence. Samorodnitsky (2007) gives a comprehensive approach to long-range dependence under heavy-tailed assumptions. We plan to address these open questions in future work.

# 8 Proof of the results

Here we prove our main results. Some proofs use technical lemmas in the Appendix. The proof is difficult because we did not find adequate a.s. convergence results adapted to the current triangular arrays setting. This led us to a systematic use of the Borel-Cantelli lemma and the bounds for the probability of the events on which the estimates do not have a satisfying behaviour. Those quantities are controlled through standard exponential inequalities or through the Markov inequality and suitable moment inequalities. From now on, we denote by $| \cdot |$ the Euclidean norm on $\mathbb{R}^2$: $|(x, y)| := (x^2 + y^2)^{1/2}$. We also denote the corresponding operator norm on the set of $2 \times 2$ matrices by $| \cdot |$, namely, $|A| := \sup\{|Az|/|z| \mid z \in \mathbb{R}^2\}$. Finally, if $X$ is a random variable, we denote by $\|X\|_p = \mathbb{E}^{1/p} [|X|^p]$ the $L_p$-norm of $X$, for any $p \ge 1$.

## 8.1 Proof of Theorem 1.

The proof will be based on the technical Lemma 1 in the Appendix. Lemma 1 controls the probability of getting small values for empirical means and empirical variances. We then use the constants depending on $\epsilon$ and $L$ defined in Lemma 1.

*Proof.* Set $\gamma \in (0, L)$. Define the following event

$$E_n(\epsilon) := \bigcap_{t=1}^{T} \{\hat{m}_t > \epsilon, \hat{\sigma}_t^2 > \epsilon\}.$$

This sequence of events depending on $n$ takes care of the properties of the estimates when they take small values. Using the union bound and Lemma 1, we have

$$\mathbb{P}(E_n(\epsilon)^c) \le 2T \exp(-nL).$$

This shows that $\sum_{n \ge 2} \mathbb{P}(E_n(\epsilon)^c) < \infty$.
We now study the consistency of $\hat{\theta}$. We only need to show that $\hat{D}_T - D_T$ and $\hat{N}_T - N_T$ go to 0 in probability or a.s. Indeed, if we use the decomposition

$$ \hat{\theta} - \theta = \hat{D}_T^{-1} (\hat{N}_T - N_T + (D_T - \hat{D}_T)D_T^{-1}N_T), $$

Assumptions **A2** and **A4** guarantee that $N_T = O(1)$ and $D_T^{-1} = O(1)$. There then also exists $\delta > 0$ sufficiently small such that $\sup_{|H|\le\delta} |(D_T + H)^{-1} - D_T^{-1}| = O(1)$.

We will then conclude that $\hat{D}_T^{-1} - D_T^{-1}$ converges to 0 in probability or a.s. The weak or strong consistency of $\hat{\theta}$ will then follow.

Now, for some $\eta > 0$, we have

$$ \mathbb{P}(|\hat{\theta} - \theta| > \eta) \le \mathbb{P}(\{|\hat{\theta} - \theta| > \eta\} \cap E_n(\epsilon)) + \mathbb{P}(E_n(\epsilon)^c). \quad (11) $$

We will use the following inequality extensively. If $a, b \ge \epsilon$ and $\kappa \in (0, 1)$, we have

$$ |\log a - \log b| = \frac{1}{\kappa} |\log a^{\kappa} - \log b^{\kappa}| \le \frac{1}{\kappa \epsilon^{\kappa}} |a - b|^{\kappa}. $$

This inequality turns the log-problems across the origin into easier moment bounds. To get this inequality, one can apply the mean value theorem, then use the inequality $|a^{\kappa} - b^{\kappa}| \le |a - b|^{\kappa}$ for $\kappa \in (0, 1)$.

On $E_n(\epsilon)$, we get

$$ |\log \hat{\sigma}_t^2 - \log \sigma_t^2| \le \frac{1}{\kappa \epsilon^{\kappa}} |\hat{\sigma}_t^2 - \sigma_t^2|^{\kappa}. $$

From Assumption **A2** and standard moment inequalities for partial sums (Lemma 5), we have

$$ \sup_{t \ge 1} \mathbb{E} |\hat{\sigma}_t^2 - \sigma_t^2|^\kappa \le \sup_{t \ge 1} \mathbb{E}^{2\kappa/k} |\hat{\sigma}_t^2 - \sigma_t^2|^{k/2} = o(1). $$

By choosing $\kappa$ small enough in the previous bounds, for any positive integer $i$, we get

$$ \sup_{t \ge 1} \mathbb{E} \left[ \left| \log \hat{\sigma}_t^2 - \log \sigma_t^2 \right|^i \mathbb{1}_{E_n(\epsilon)} \right] = o(1) \quad (12) $$

and similarly

$$ \sup_{t \ge 1} \mathbb{E} \left[ \left| \log \hat{\mu}_t - \log \mu_t \right|^i \mathbb{1}_{E_n(\epsilon)} \right] = o(1). \quad (13) $$

On the event $E_n(\epsilon)$, $T^{-1} \sum_{t=1}^T [\log \hat{\sigma}_t^2 - \log \sigma_t^2] \to 0$ in $\mathbb{L}^1$ and then in probability when $n \to \infty$. Moreover, setting $C := \max(\sup_{t \ge 1} |\log \sigma_t^2|, \sup_{t \ge 1} |\log \mu_t|)$, we have

$$ \begin{aligned} & |\log \hat{\mu}_t \log \hat{\sigma}_t^2 - \log \mu_t \log \sigma_t^2| \\ & \le |\log \hat{\mu}_t - \log \mu_t| \cdot |\log \hat{\sigma}_t^2 - \log \sigma_t^2| + |\log \sigma_t^2| \cdot |\log \hat{\mu}_t - \log \mu_t| + |\log \mu_t| \cdot |\log \hat{\sigma}_t^2 - \log \sigma_t^2| \\ & \le \frac{1}{2} |\log \hat{\sigma}_t^2 - \log \sigma_t^2|^2 + \frac{1}{2} |\log \hat{\mu}_t - \log \mu_t|^2 + C |\log \hat{\mu}_t - \log \mu_t| + C |\log \hat{\sigma}_t^2 - \log \sigma_t^2|. \end{aligned} $$
From (12) and (13), we also deduce that

$$ \frac{1}{T} \sum_{t=1}^{T} \mathbb{E}1_{E_n(\epsilon)} |\log \hat{\mu}_t \log \hat{\sigma}_t^2 - \log \mu_t \log \sigma_t^2| = o(1). $$

Hence $\hat{N}_T - N_T \to 0$ in probability on the event $E_n(\epsilon)$. Similarly, $\hat{D}_T - D_T \to 0$ in probability on the event $E_n(\epsilon)$. This leads to the weak consistency of the regression estimates.

For any $\eta > 0$,

$$ \sum_{n \ge 1} \mathbb{P} \left( \left\{ \left| \hat{\theta} - \theta \right| > \eta \right\} \cap E_n(\epsilon)^c \right) \le \sum_{n \ge 1} \mathbb{P} (E_n(\epsilon)^c) < \infty $$

and the Borel-Cantelli lemma ensures that $|\hat{\theta} - \theta| \cdot 1_{E_n(\epsilon)^c}$ converges to 0 a.s. To prove the almost sure convergence, it is enough to prove that $\hat{N}_T - N_T \to 0$ and $\hat{D}_T - D_T \to 0$ a.s. on the event $E_n(\epsilon)$. Using the same strategy as for the convergence in probability, it suffices to show that

$$ \lim_{T \to \infty} \frac{1}{T} \sum_{t=1}^{T} \{ |\hat{\sigma}_t^2 - \sigma_t^2| + |\hat{\mu}_t - \mu_t| \} = 0 \quad \text{a.s.} \qquad (14) $$

Because the function $x \mapsto |x|^{k/2}$ (where $k > 2$ is given in Assumption A2 is convex, we get

$$ M_T := \left( \frac{1}{T} \sum_{t=1}^{T} |\hat{\sigma}_t^2 - \sigma_t^2| \right)^{k/2} \le \frac{1}{T} \sum_{t=1}^{T} |\hat{\sigma}_t^2 - \sigma_t^2|^{k/2}. $$

From moment inequalities for partial sums of independent observations (Lemma 5), we get $EM_T = O(n^{-k/4})$ and by assumption, if $k > 4$, we get $\sum_{n \ge 1} EM_T < \infty$. Setting $P_T = |\hat{\mu}_t - \mu_t|$, similar arguments yield $\sum_{n \ge 1} EP_T < \infty$. The Borel-Cantelli lemma ensures that both $M_T$ and $P_T$ go to zero a.s. as $n \to \infty$ and then the validity of (14), which concludes the proof. $\square$

**Proof of Theorem 2.** The proof is similar to that of Theorem 1. The main difference is in proving the convergence of $\hat{D}_T - D_T$ or $\hat{N}_T - N_T$ to 0.

*Proof.* Let $A_t = \{\hat{\mu}_t > \epsilon, \hat{\sigma}_t^2 > \epsilon\}$ where $\epsilon < 1$ satisfies the condition in Lemma 1.2. Using arguments similar to those in the proof of Theorem 1, one can show that

$$ \frac{1}{T} \sum_{t=1}^{T} 1_{A_t} \left[ |\log \hat{\sigma}_t^2 - \log \sigma_t^2| + |\log \hat{\mu}_t - \log \mu_t| + |\log^2 \hat{\mu}_t - \log^2 \mu_t| + |\log \hat{\mu}_t \log \hat{\sigma}_t^2 - \log \mu_t \log \sigma_t^2| \right] $$

converges to 0 either in probability or a.s., depending on $k$. From Lemma 1,

$$ \sum_{n \ge 1} \frac{1}{T} \sum_{t=1}^{T} \mathbb{P}(A_t^c) < \infty $$
and from the Borel-Cantelli lemma, we deduce that $\lim_{n \to \infty} \frac{1}{T} \sum_{t=1}^{T} 1_{A_t^c} = 0$ a.s. Moreover, from Lemma 1, the sequences $(\mu_t)_{t \ge 1}$ and $(\sigma_t^2)_{t \ge 1}$ are upper- and lower-bounded. To show the required convergence, it suffices to show that the following quantities converge a.s. to 0:

$$ \frac{1}{T} \sum_{t=1}^{T} 1_{0 < \hat{\sigma}_t^2 < \epsilon} \cdot |\log \hat{\sigma}_t^2|, $$

$$ \frac{1}{T} \sum_{t=1}^{T} 1_{0 < \hat{\sigma}_t^2 < \epsilon} \cdot |\log \hat{\sigma}_t^2 \cdot \log \hat{\mu}_t|, $$

$$ \frac{1}{T} \sum_{t=1}^{T} 1_{0 < \hat{\mu}_t < \epsilon} \cdot |\log \hat{\sigma}_t^2 \cdot \log \hat{\mu}_t|, $$

$$ \frac{1}{T} \sum_{t=1}^{T} 1_{0 < \hat{\mu}_t < \epsilon} \cdot |\log \hat{\mu}_t|, $$

$$ \frac{1}{T} \sum_{t=1}^{T} 1_{0 < \hat{\mu}_t < \epsilon} \cdot |\log^2 \hat{\mu}_t|. $$

We recall that $\log \hat{\mu}_t$ and $\log \hat{\sigma}_t^2$ are implicitly defined respectively by $\log \hat{\mu}_t \cdot 1_{\hat{\mu}_t > 0}$ and $\log \hat{\sigma}_t^2 \cdot 1_{\hat{\sigma}_t^2 > 0}$. First, because $\log(\epsilon) < 0$, we have that

$$ \log \hat{\mu}_t \cdot 1_{0 < \hat{\mu}_t < \epsilon} = \log \hat{\mu}_t \cdot 1_{0 < \hat{\mu}_t < \epsilon} \cdot 1_{\cup_{j=1}^n \{X_{j,t} > 0\}} \ge \sum_{j=1}^n \log (X_{j,t}/n) \cdot 1_{X_{j,t} > 0} \cdot 1_{0 < \hat{\mu}_t < \epsilon}. $$

We then get the bound

$$ |\log \hat{\mu}_t| \cdot 1_{\hat{\mu}_t > 0} \le |\log \hat{\mu}_t| \cdot 1_{\hat{\mu}_t > \epsilon} - \sum_{j=1}^{n} \log(X_{j,t}/n) \cdot 1_{X_{j,t} > 0} \cdot 1_{0 < \hat{\mu}_t < \epsilon}. \quad (15) $$

Using the same type of arguments for

$$ \hat{\sigma}_t^2 = \frac{1}{2n^2} \sum_{1 \le i \ne j \le n} |X_{i,t} - X_{j,t}|^2, $$

we get the bound

$$ |\log \hat{\sigma}_t^2| \cdot 1_{\hat{\sigma}_t^2 > 0} \le |\log \hat{\sigma}_t^2| \cdot 1_{\hat{\sigma}_t^2 > \epsilon} - \sum_{1 \le i \ne j \le n} \log \left( \frac{|X_{i,t} - X_{j,t}|^2}{2n^2} \right) \cdot 1_{X_{i,t} \ne X_{j,t}} \cdot 1_{0 < \hat{\sigma}_t^2 < \epsilon}. \quad (16) $$

We are going to show that

$$ Z_T := \frac{1}{T} \sum_{t=1}^{T} 1_{0 < \hat{\sigma}_t^2 < \epsilon} \cdot \log \hat{\sigma}_t^2 \cdot \log \hat{\mu}_t \to_{n \to \infty} 0 \quad (17) $$

in probability or a.s. The idea for all the other quantities is exactly the same and details will be omitted.

Let $p = q/2 > 1$ by Assumption A5 and $\alpha = 1/(2p)$. From (15), there exists $C_q > 0$ such that

$$ \| \log \hat{\mu}_t \cdot 1_{\hat{\mu}_t > 0} \|_{2p} \le C_q \cdot \mathbb{E}^{\alpha}(\hat{\mu}_t) + \sum_{j=1}^{n} \| \log \frac{X_{j,t}}{n} 1_{X_{j,t} > 0} \|_{2p}. $$

From A5, we deduce the existence of a constant $D > 0$ such that the behaviour of $\hat{\mu}_t$ around the origin is not too explosive:

$$ \| \log \hat{\mu}_t \cdot 1_{\hat{\mu}_t > 0} \|_{2p} \le D (1 + n \log(n)). \quad (18) $$
From (16) and the Hölder inequality, we deduce that

$$ \begin{aligned} & \| \log \hat{\mu}_t \cdot \log \hat{\sigma}_t^2 \cdot \mathbb{1}_{0<\hat{\sigma}_t^2<\epsilon} \|_1 \\ & \le \sum_{1 \le i \ne j \le n} \mathbb{P} (\hat{\sigma}_t^2 < \epsilon)^{1-1/p} \cdot \| \log \hat{\mu}_t \cdot \mathbb{1}_{\hat{\mu}_t > 0} \|_{2p} \cdot \| \log \left( \frac{|X_{i,t} - X_{j,t}|^2}{2n^2} \right) \cdot \mathbb{1}_{X_{i,t} \ne X_{j,t}} \|_{2p}. \end{aligned} $$

From Assumption A5, inequality (18) and Lemma 1, it is not very difficult to conclude that

$$ \sum_{n \ge 1} \frac{1}{T} \sum_{t=1}^{T} \| \log \hat{\mu}_t \cdot \log \hat{\sigma}_t^2 \cdot \mathbb{1}_{0<\hat{\sigma}_t^2<\epsilon} \|_1 < \infty. $$

From the Borel-Cantelli lemma, we get (17). Remembering the discussion at the beginning of the proof, one can deduce the consistency properties of $\hat{\theta}$. □

## 8.2 Proof of Theorem 3.

*Proof.* The second derivatives of the function $\phi_{\theta} = (\phi_{1,\theta}, \phi_{2,\theta})$ are:

$$ \begin{aligned} \frac{\partial^2 \phi_{1,\theta}}{\partial x^2}(x,y) &= -\frac{2}{y-x^2} - \frac{4x^2}{(y-x^2)^2} + \frac{\theta_2}{x^2}, \\ \frac{\partial^2 \phi_{1,\theta}}{\partial y^2}(x,y) &= \frac{-1}{(y-x^2)^2}, \quad \frac{\partial^2 \phi_{1,\theta}}{\partial x \partial y}(x,y) = \frac{2x}{(y-x^2)^2}, \\ \frac{\partial^2 \phi_{2,\theta}}{\partial x^2}(x,y) &= \frac{-4}{y-x^2} - \frac{\log(y-x^2)}{x^2} - \frac{2\log x}{y-x^2} - \frac{4x^2 \log(x)}{(y-x^2)^2} + \frac{\theta_1 - 2\theta_2}{x^2} + \frac{2\theta_2 \log(x)}{x^2}, \\ \frac{\partial^2 \phi_{2,\theta}}{\partial y^2}(x,y) &= -\frac{\log(x)}{(y-x^2)^2}, \quad \frac{\partial^2 \phi_{2,\theta}}{\partial x \partial y}(x,y) = \frac{1}{x(y-x^2)} + \frac{2x \log(x)}{(y-x^2)^2}. \end{aligned} $$

To prove Theorem 3.1, we first show that

$$ \sqrt{nT}(U_T - V_T) = o_P(1). \quad (19) $$

Using a Taylor-Lagrange expansion of order 2, we have

$$ \sqrt{nT}(U_T - V_T) = \left( \frac{1}{2} \frac{\sqrt{n}}{\sqrt{T}} \sum_{t=1}^{T} (\hat{m}_t - m_t)' H_{1,t} (\hat{m}_t - m_t), \frac{1}{2} \frac{\sqrt{n}}{\sqrt{T}} \sum_{t=1}^{T} (\hat{m}_t - m_t)' H_{2,t} (\hat{m}_t - m_t) \right), (20) $$

where

$$ H_{1,t} = \nabla^{(2)} \phi_{1,\theta}(\tilde{m}_t), \quad H_{2,t} = \nabla^{(2)} \phi_{2,\theta}(\tilde{m}_t), $$

for some random points $\tilde{m}_t, \tilde{m}_t$ in the segment $[m_t, \hat{m}_t]$. It is only necessary to show (19) on the event

$$ E(\epsilon) = \cap_{t=1}^{T} \{ \hat{\mu}_t > \epsilon, \hat{\sigma}_t^2 > \epsilon \}, $$

which has been defined in the proof of Theorem 1 and which satisfies $\mathbb{P}(E(\epsilon)^c) \le 2T \exp(-nL)$ for some $L > 0$. From Lemma 1.1, on $E(\epsilon)$, we have, setting $w = \min(\epsilon, s)$,

$$ \min_{1 \le t \le T} \min (\mu_t, m_{2,t}, \sigma_t^2, \hat{\mu}_t, \hat{m}_{2,t}, \hat{\sigma}_t^2) \ge w. $$
Moreover, on the same event $E(\epsilon)$, using convexity of the square function, we have for $1 \le t \le T$,

$$ \min (\tilde{m}_{2,t} - \tilde{m}_{1,t}^2, \tilde{m}_{2,t} - \tilde{m}_{1,t}^2) \ge \min (\hat{\sigma}_t^2, \sigma_t^2) \ge w. $$

Also $\min (\tilde{m}_{1,t}, \tilde{m}_{1,t}) \ge w$. Moreover, there exists a positive constant $\tilde{C}$ such that for $1 \le t \le T$ and $i = 1, 2$,

$$ \max (\tilde{m}_{i,t}, \tilde{m}_{i,t}) \le \tilde{C} (1 + |\hat{m}_{i,t} - m_{i,t}|). $$

It is then sufficient to show that for $i, i' \in \{1, 2\}$,

$$ \frac{\sqrt{n}}{\sqrt{T}} \sum_{t=1}^{T} h_{i,i',t} |\hat{m}_{i,t} - m_{i,t}| \cdot |\hat{m}_{i',t} - m_{i',t}| = o_{\mathbb{P}}(1), \qquad (21) $$

with

$$ h_{i,i',t} = \max_{\ell=1,2} \sup_{\substack{w \le x \le \tilde{C}(1+|\hat{\mu}_t - \mu_t|) \\ w \le z \le \tilde{C}(1+|\hat{m}_{2,t} - m_{2,t}|)}} \left| \frac{\partial^2 \phi_{\ell,\theta}}{\partial x^i \partial y^{i'}} (x, z + x^2) \right|. $$

A careful inspection of the second partial derivatives of $\phi_\theta$ shows that for any $\delta > 0$, there exists a constant $C'_\delta > 0$ such that

$$ h_{1,1,t} \le C'_\delta [1 + |\hat{\mu}_t - \mu_t|^{2+\delta} + |\hat{m}_{2,t} - m_{2,t}|^\delta], $$

$$ h_{2,2,t} \le C'_{\delta} [1 + |\hat{\mu}_t - \mu_t|^\delta], \quad h_{1,2,t} \le C'_{\delta} [1 + |\hat{\mu}_t - \mu_t|^{1+\delta}]. $$

Take $\delta \in (0, 1)$ such that $4 + \delta < k$. Using the moment bound given by Lemma 5, it is not difficult to get that

$$ \frac{\sqrt{n}}{\sqrt{T}} \sum_{t=1}^{T} h_{i,i',t} |\hat{m}_{i,t} - m_{i,t}| \cdot |\hat{m}_{i',t} - m_{i',t}| = O_{\mathbb{P}}(\sqrt{T/n}) $$

which yields (21) and then (19).

It now remains to show that

$$ W_T := \sqrt{nT} C_T V_T \Rightarrow \mathcal{N}_2(0, I_2). \qquad (22) $$

We recall that $I_2$ denotes the $2 \times 2$ identity matrix.

This will follow from the Lindeberg-Feller central limit theorem for triangular arrays (Billingsley, 2013, Theorem 18.1). Since $\mathrm{Var}(W_T) = I_2$, we only have to check the Lindeberg condition. Setting $Z_{n,j,t} = \frac{1}{\sqrt{nT}} C_T J_{\phi_\theta}(m_t)(Y_{j,t} - m_t)$ and using our assumptions, we have $C_T = O(1)$ and $\max_{1 \le t \le T} |J_{\phi_\theta}(m_t)| = O(1)$. For $\delta > 0$ such that $2(2 + \delta) \le k$, we get $\sum_{j=1}^n \sum_{t=1}^T \mathbb{E}|Z_{n,j,t}|^{2+\delta} = O((nT)^{-\delta/2}) = o(1)$. This proves (22).

To prove Theorem 3.2, we again use the decomposition (20) and we first show that for $\ell = 1, 2$,

$$ \frac{\sqrt{n}}{\sqrt{T}} \sum_{t=1}^{T} (\hat{m}_t - m_t)' [H_{\ell,t} - K_{\ell,t}] (\hat{m}_t - m_t) = o_{\mathbb{P}}(1), \qquad (23) $$
where for $\ell = 1, 2$ and $t \ge 1$, $K_{\ell,t} = \nabla^{(2)}\phi_{\ell,\theta}(m_t)$. It is only necessary to show (23) on the event $E_n(\epsilon)$. To this end, we will use the bounds

$$ \max \{|\tilde{\mu}_t - \mu_t|, |\tilde{\mu}_t - \mu_t|\} \le |\hat{\mu}_t - \mu_t|, $$

and

$$ \max \{|\tilde{\sigma}_t^2 - \sigma_t^2|, |\tilde{\sigma}_t^2 - \sigma_t^2|\} \le |\hat{m}_{2,t} - m_{2,t}| + \left(2 \sup_{s \ge 1} \mu_s + |\hat{\mu}_t - \mu_t|\right) |\hat{\mu}_t - \mu_t|, $$

where we set $\tilde{\sigma}_t^2 := \tilde{m}_{2,t} - \tilde{m}_{1,t}^2$ and $\tilde{\sigma}_t^2 := \tilde{m}_{2,t} - \tilde{m}_{1,t}^2$. We will also use the fact that if $h : (\epsilon, \infty) \to \mathbb{R}_+$ is any mapping of the form $h(z) = 1/z$, $h(z) = 1/z^2$ or $h(z) = \log(z)$, then for any $\kappa \in (0, 1)$, there exists $C_\kappa > 0$ such that

$$ |h(z) - h(z')| \le C_\kappa |z - z'|^\kappa. $$

Assertion (23) can be obtained from these facts and a careful analysis of the second derivatives of $\phi_\theta$. We will illustrate this by showing that

$$ \frac{\sqrt{n}}{\sqrt{T}} \sum_{t=1}^{T} (\hat{m}_{2,t} - m_{2,t})^2 \left[ \frac{\partial^2 \phi_{1,\theta}}{\partial y^2}(\tilde{m}_t) - \frac{\partial^2 \phi_{1,\theta}}{\partial y^2}(m_t) \right] = o_P(1). \quad (24) $$

For $\kappa \in (0, 1)$ small enough to satisfy $4 + 2\kappa \le k$, there exists $D_\kappa > 0$ such that

$$ \left| \frac{\tilde{\mu}_t}{\tilde{\sigma}_t^4} - \frac{\mu_t}{\sigma_t^4} \right| \le C \left[ |\tilde{\mu}_t - \mu_t| + |\tilde{\sigma}_t^2 - \sigma_t^2|^\kappa \right]. $$

To show (24), it is sufficient to show that

$$ \mathbb{E} \left[ (\hat{m}_{2,t} - m_{2,t})^2 |\hat{m}_{i,t} - m_{i,t}|^{\kappa_i} \right] = O \left( n^{-1-\frac{\kappa_i}{2}} \right), $$

when $(i, \kappa_i) = (2, \kappa)$, $(1, \kappa)$ or $(1, 2\kappa)$. But this point easily follows from Lemma 5 and Hölder's inequality. Assertions similar to (24) can be obtained in the same way for the other second-order partial derivatives of $\phi_\theta$. Details are omitted. Then

$$ \sqrt{nT}(U_T - V_T) = \left( \frac{\sqrt{n}}{\sqrt{T}} \sum_{t=1}^{T} (\hat{m}_t - m_t)' K_{1,t} (\hat{m}_t - m_t), \frac{\sqrt{n}}{\sqrt{T}} \sum_{t=1}^{T} (\hat{m}_t - m_t)' K_{2,t} (\hat{m}_t - m_t) \right) + o_P(1). $$

One applies Lemma 2 to conclude the proof. $\square$

## 8.3 Proof of Corollary 2.

The result will be consequence of Theorem 3 and of Lemmas 7 and 4 in the Appendix.

*Proof.* Set $\hat{A}_n := \hat{C}_T \hat{D}_T$ and $A_n := C_T \hat{D}_T$. Using Lemma 7 in the Appendix, we have $\hat{C}_T - C_T = o_P(1)$ and from the proof of Theorem 1, we know that $\hat{D}_T - D_T = o_P(1)$ and both $D_T$ and $D_T^{-1}$ are $O(1)$. We then easily conclude that $\hat{A}_n - A_n = o_P(1)$.

Let us also define $\hat{B}_n := \sqrt{T/n}\hat{C}_T\hat{E}_T$ if $T/n = O(1)$ and $\hat{B}_n := 0$ if $T/n = o(1)$. Once again, using Lemma 7, it is easy to show that $\hat{B}_n - B_n = o_P(1)$ where $B_n := \sqrt{T/n}C_T E_T$ if $T/n = O(1)$
and $B_n := 0$ if $T/n = o(1)$. Setting $W_n = \sqrt{nT}(\hat{\theta} - \theta)$, we already know from Theorem 3 that $A_nW_n - B_n \Rightarrow N_2(0, I_2)$. Moreover, using the decomposition

$$ 
\begin{aligned}
Z_n &:= \hat{A}_n W_n - \hat{B}_n \\
&= A_n W_n - B_n \\
&+ (\hat{A}_n - A_n) A_n^{-1} (A_n W_n - B_n) + (\hat{A}_n A_n^{-1} - I_2) B_n + B_n - \hat{B}_n.
\end{aligned}
$$

We conclude from Lemma 7 that $Z_n \Rightarrow N_2(0, I_2)$. The result is now a consequence of Lemma 4. □

# 9 Appendix

This Appendix states and proves several technical lemmas used to prove our main results.

Lemma 1 below shows that Assumption A3 leads to a lower bound for the population mean and the population variance. It also shows that small values of the sample mean and sample variance occur with very small probability if $n$ is large, as noted at the beginning of the proofs section 8.

**Lemma 1.** Suppose that Assumption A3 holds, i.e., that there exists some $\delta > 0$ such that

$$ p(\delta) := \inf_{t \ge 1} \mathbb{P}(X_{1,t} > \delta) > 0 \quad \text{and} \quad q(\delta) = \inf_{t \ge 1} \mathbb{P}(|X_{1,t} - X_{2,t}| > \delta) > 0. $$

Then

1. There exists $s > 0$ such that $\inf_{t \ge 1} \mu_t \ge s$ and $\inf_{t \ge 1} \sigma_t^2 \ge s$.

2. Let $0 < \epsilon < \min(\delta p(\delta), \delta^2 q(\delta)/4)$. Then there exists $L > 0$ such that

$$ \sup_{t \ge 1} \max \{ \mathbb{P}(\hat{\mu}_t \le \epsilon), \mathbb{P}(\hat{\sigma}_t^2 \le \epsilon) \} \le \exp(-Ln). $$

*Proof.* To prove point 1, we use the inequality $x^i \ge \delta^i \mathbb{1}_{\{x \ge \delta\}}$ for any $x \ge 0$ and $i = 1, 2$, which leads to the lower bounds

$$ \mu_t \ge \delta \mathbb{P}(X_{1,t} > \delta) \ge \delta p(\delta), $$

$$ \sigma_t^2 = \frac{1}{2} \mathbb{E}(|X_{1,t} - X_{2,t}|^2) \ge \frac{\delta^2}{2} \mathbb{P}(|X_{1,t} - X_{2,t}| \ge \delta) \ge \frac{\delta^2}{2} q(\delta). $$

To prove point 2, we observe that

$$ \hat{\mu}_t \ge \frac{\delta}{n} \sum_{j=1}^{n} \mathbb{1}_{X_{j,t} > \delta}. $$

Then, setting $Z_{j,t} := \mathbb{1}_{X_{j,t} > \delta} - \mathbb{P}(X_{j,t} > \delta)$,

$$ \mathbb{P}(\hat{\mu}_t \le \epsilon) \le \mathbb{P}\left(\frac{1}{n} \sum_{j=1}^{n} Z_{j,t} \le \frac{\epsilon}{\delta} - p(\delta)\right). $$

Using Hoeffding's exponential inequality (Hoeffding, 1963, inequality (2.6)), we conclude that $\sup_{t \ge 1} \mathbb{P}(\hat{\mu}_t \le \epsilon) \le \exp(-nL')$ for some $L' > 0$. $L'$ can be taken here as $(\epsilon/\delta - p(\delta))^2$. To
control the probability of small values of the empirical variance, we use an immediate consequence of Lagrange's identity:

$$ \hat{\sigma}_t^2 = \frac{1}{2n^2} \sum_{1 \le i \ne j \le n} (X_{i,t} - X_{j,t})^2 . $$

This leads to the lower bound

$$ \hat{\sigma}_t^2 \ge \frac{\delta^2}{2n^2} \sum_{1 \le i \ne j \le n} 1_{|X_{i,t} - X_{j,t}| > \delta} . $$

Next, set $Z_{i,j}^{(t)} := 1_{|X_{i,t} - X_{j,t}| > \delta} - \mathbb{P}(|X_{i,t} - X_{j,t}| > \delta)$. Then

$$ \mathbb{P}(\hat{\sigma}_t^2 \le \epsilon) \le \mathbb{P} \left( \frac{1}{n(n-1)} \sum_{1 \le i \ne j \le n} Z_{i,j}^{(t)} \le \frac{2n^2\epsilon}{\delta^2 n(n-1)} - q(\delta) \right) . $$

Since $2n^2/(n(n-1)) \le 4$ for every $n \ge 2$, an application of Hoeffding's inequality for U-statistics (Hoeffding, 1963, bound (5.7)) leads to

$$ \mathbb{P}(\hat{\sigma}_t^2 \le \epsilon) \le \exp(-2[n/2]t^2) $$

where $[n/2]$ denotes the integer-part of $n/2$ and $t = q(\delta) - \frac{4\epsilon}{\delta^2}$. It is then always possible to get the conclusion of the Lemma for a suitably chosen constant $L \in (0, L')$. $\square$

**Lemma 2.** Suppose that Assumption A1 holds true and for any $j \ge 1$, the variables $X_{j,1}, \dots, X_{j,T}$ are mutually independent. Let $\{J_t : t \in \mathbb{N}\}$ be a family of $2 \times 2$ matrices such that $\sup_{t \ge 1} |J_t| < \infty$. Suppose that $\sup_{t \ge 1} \mathbb{E}X_{1,t}^{4(1+\kappa)} < \infty$ for some $\kappa \in (0, 1)$. Then

$$ \frac{1}{T} \sum_{t=1}^{T} (\hat{m}_t - m_t)' J_t (\hat{m}_t - m_t) = \frac{1}{nT} \sum_{t=1}^{T} \mathrm{Tr}(J_t \mathrm{Var}(Y_{1,t})) + O_P \left( \left( \frac{T}{n} \right)^{\frac{1}{2} - \frac{\kappa}{1+\kappa}} \frac{1}{n^{\frac{2\kappa}{1+\kappa}} \sqrt{nT}} + \frac{1}{\sqrt{n} \sqrt{nT}} \right), $$

where $Y_{j,t} = (X_{j,t}, X_{j,t}^2)'$ for $1 \le j \le n$ and $1 \le t \le T$.

*Proof.* Setting $\bar{Y}_{j,t} := Y_{j,t} - \mathbb{E}Y_{j,t}$ and

$$ A := \frac{1}{n^2 T} \sum_{t=1}^{T} \sum_{j=1}^{n} \bar{Y}_{j,t}' J_t \bar{Y}_{j,t}, $$

$$ B := \frac{2}{n^2 T} \sum_{t=1}^{T} \sum_{1 \le i < j \le n} \bar{Y}_{i,t}' J_t \bar{Y}_{j,t}, $$

we have

$$ \frac{1}{T} \sum_{t=1}^{T} (\hat{m}_t - m_t)' J_t (\hat{m}_t - m_t) = A + B. $$

It is straightforward to get

$$ \mathbb{E}(A) = \frac{1}{nT} \sum_{t=1}^{T} \mathrm{Tr}(J_t \mathrm{Var}(Y_{1,t})). $$
Moreover, setting $p = 1 + \kappa$, one can apply Lemma 5 to get

$$ \|A - \mathbb{E}(A)\|_p = O\left(n^{-1-\frac{\kappa}{1+\kappa}}T^{-\frac{\kappa}{1+\kappa}}\right). $$

Then

$$ \begin{aligned} \mathbb{E}(B^2) &= \frac{4}{n^4 T^2} \sum_{t=1}^{T} \sum_{1 \le i < j \le n} \mathbb{E}(|\bar{Y}'_{i,t} J_t \bar{Y}_{j,t}|^2) \\ &\le \frac{2n(n-1)}{n^4 T} \sup_{t \ge 1} |J_t| \times \sup_{t \ge 1} \mathbb{E}^2(|\bar{Y}_{1,t}|^2). \end{aligned} $$

$$ \text{Hence } B = O_{\mathbb{P}}\left(\frac{1}{\sqrt{n\sqrt{nT}}}\right). \quad \square $$

**Lemma 3.** Assume that $(Z_n)_{n \in \mathbb{N}}$ is a sequence of random vectors of $\mathbb{R}^2$ such that $Z_n \Rightarrow \mathcal{N}(0, I_2)$. Then for any sequence $(u_n)_{n \in \mathbb{N}}$ of deterministic vectors of $\mathbb{R}^2$ such that $|u_n| = 1$, we have

$$ u'_{n}Z_{n} \Rightarrow \mathcal{N}(0, 1). $$

*Proof.* It is clear that the sequence $(u'_{n}Z_{n})_{n \in \mathbb{N}}$ is bounded in probability (i.e., tight). We simply have to show that any subsequence has a subsubsequence converging to the distribution $\mathcal{N}(0, 1)$. Let $(u'_{\phi_n} Z_{\phi_n})_{n \in \mathbb{N}}$ be a given subsequence. Since the unit sphere is compact, one can find a subsubsequence $(u_{\psi_n})_{n \in \mathbb{N}}$ with $\psi_n \in \{\phi_k : k \in \mathbb{N}\}$ for all $n \in \mathbb{N}$ and such that $\lim_{n \to \infty} u_{\psi_n} = u$ for some $u \in \mathbb{R}^2$ such that $|u| = 1$. Since

$$ u'_{\psi_n} Z_{\psi_n} = (u_{\psi_n} - u)' Z_{\psi_n} + u' Z_{\psi_n}, $$

our assumption guarantees that $u'Z_{\psi_n} \Rightarrow \mathcal{N}(0, 1)$ and we also know that $(u_{\psi_n} - u)' Z_{\psi_n}$ will converge in probability to 0 as $n \to \infty$. From Slutsky's lemma, we conclude that $u'_{\psi_n} Z_{\psi_n} \Rightarrow \mathcal{N}(0, 1)$. $\square$

**Lemma 4.** Suppose $Z_n = \hat{A}_n W_n - \hat{B}_n \Rightarrow \mathcal{N}_2(0, I_2)$ with $\hat{A}_n$ a random real matrix of size $2 \times 2$, $W_n$ a random real vector of $\mathbb{R}^2$ and $\hat{B}_n$ a random real vector of $\mathbb{R}^2$. Suppose furthermore that $\hat{A}_n - A_n = o_{\mathbb{P}}(1)$, with $A_n$ deterministic and such that $A_n = O(1)$ and $A_n^{-1} = O(1)$. Setting $H_n := \hat{A}_n^{-1}(\hat{A}_n^{-1})'$ and $m_n = \hat{A}_n^{-1}\hat{B}_n$, we have for $i = 1, 2$,

$$ \frac{W_n(i) - m_n(i)}{\sqrt{H_n(i, i)}} \Rightarrow \mathcal{N}(0, 1). $$

*Proof.* To show first that

$$ \hat{A}_n^{-1} - A_n^{-1} = o_{\mathbb{P}}(1), \quad (25) $$

we write

$$ \hat{A}_n^{-1} = (A_n - (A_n - \hat{A}_n))^{-1} = \sum_{k \ge 0} (A_n^{-1}(A_n - \hat{A}_n))^k A_n^{-1}, $$

which is bounded in probability because $\hat{A}_n - A_n \to 0$ in probability. We then deduce that

$$ \hat{A}_n^{-1} - A_n^{-1} = \hat{A}_n^{-1} (A_n - \hat{A}_n) A_n^{-1} = o_P(1), $$
which proves (25). Next, for $i = 1, 2$, we set $\hat{v}_{n,i} := (\hat{A}_n^{-1})'e_i$, where $e_1 := (1,0)'$, $e_2 := (0,1)'$. Define also $v_{n,i} := (A_n^{-1})'e_i$, $i = 1, 2$. Since $e_i = A_n'v_{n,i}$ and $A_n = O(1)$, $v_{n,i}$ cannot approximate 0. Moreover, $\hat{v}_{n,i} - v_{n,i} = o_P(1)$. Using also that $A_n^{-1} = O(1)$, $v_{n,i}$ is bounded and we deduce that

$$ \frac{\hat{v}_{n,i}}{|\hat{v}_{n,i}|} - \frac{v_{n,i}}{|v_{n,i}|} = o_P(1). $$

We then conclude that

$$ \frac{\hat{v}'_{n,i} Z_n}{|\hat{v}_{n,i}|} - \frac{v'_{n,i} Z_n}{|v_{n,i}|} = o_P(1). $$

Using Lemma 3, we know that $v_{n,i}'Z_n/|v_{n,i}| \Rightarrow \mathcal{N}(0,1)$. Since

$$ \frac{\hat{v}'_{n,i} Z_n}{|\hat{v}_{n,i}|} = \frac{W_n(i) - m_n(i)}{|\hat{v}_{n,i}|}, $$

we only have to show that $|v_{n,i}| = \sqrt{H_n(i, i)}$ to finish the proof. □

The following lemma is useful to bound moments of partial sums. It is a very special case of the Burkholder inequality for martingales (Hall and Heyde, 2014, Theorem 2.10).

**Lemma 5.** If $\xi_1, \dots, \xi_n$ are independent random variables with 0 mean and a finite moment of order $p > 1$, then there exists a constant $C_p > 0$ such that

$$ \left\| \sum_{i=1}^{n} \xi_i \right\|_p \le C_p \times \begin{cases} n^{1/p} \max_{1 \le i \le n} \|\xi_i\|_p & \text{if } p \in (1, 2), \\ n^{1/2} \max_{1 \le i \le n} \|\xi_i\|_{p/2}^{1/2} & \text{if } p \ge 2. \end{cases} $$

**Corollary 3.** Suppose that Assumption A1 is fulfilled and that for some $k > 4$, $\sup_{t \ge 1} \mathbb{E}|X_{1,t}|^k < \infty$. Define for $i = 1, 2, 3, 4$, $S_{i,t} := \frac{1}{n} \sum_{j=1}^n (X_{j,t}^i - \mathbb{E}[X_{j,t}^i])$. Suppose that for some integer $d \ge 1$ and $i = 1, \dots, 4$, $\kappa_i$ is a positive real number, $\alpha_i \in \{1, \dots, 4\}$ with $\sum_{i=1}^d \alpha_i \kappa_i \le k$. Then

$$ \sup_{t \ge 1} \mathbb{E} \left[ \prod_{i=1}^{d} |S_{\alpha_i, t}|^{\kappa_i} \right] = O \left( n^{- \sum_{i=1}^{d} \kappa_i \min\{1/2, 1-\alpha_i/k\}} \right). $$

*Proof.* Setting $p_i := k/(\alpha_i \kappa_i)$ for $1 \le i \le d$, we have $\sum_{i=1}^d 1/p_i \le 1$ and applying the Hölder inequality, we get

$$ \mathbb{E} \left[ \prod_{i=1}^{d} |S_{\alpha_i, t}|^{\kappa_i} \right] \le \prod_{i=1}^{d} \mathbb{E}^{\frac{1}{p_i}} \left[ |S_{\alpha_i, t}|^{\frac{k}{\alpha_i}} \right]. $$

The result then follows from Lemma 5. □

**Lemma 6.** Given the assumptions of Theorem 3, let $\Gamma_T$ be the symmetric, non-negative definite matrix

$$ \Gamma_T := \frac{1}{T} \sum_{t=1}^{T} J_{\phi_{\theta}}(m_t) \mathrm{Var}(Y_{1,t}) J_{\phi_{\theta}}(m_t)' $$

Let $\lambda_T$ be the smallest eigenvalue of $\Gamma_T$ and

$$ A_T = \frac{1}{T} \sum_{t=1}^{T} J_{\phi_{\theta}}(m_t). $$
1. There exists $\xi > 0$ such that

$$ \lambda_T \geq \xi \min_{1 \leq t \leq T} \det[\mathrm{Var}(Y_{1,t})] \det^2(A_T). $$

2. Suppose that $X_{1,t}$ follows a Poisson distribution with parameter $\mu_t = \mu(t/T)$ with $\mu: [0, 1] \to (0, \infty)$ a continuous function. Then there exists $\zeta > 0$ such that for $T$ large enough, $\lambda_T \geq \zeta$.

*Proof.* 1. Let $\gamma_t$ be the smallest eigenvalue of $\mathrm{Var}(Y_{1,t})$ and $\beta_T$ the smallest eigenvalue of the matrix $A_T A'_T$. If $x \in \mathbb{R}^2$, then

$$ \begin{aligned} x'\Gamma_T x &\geq \frac{1}{T} \sum_{t=1}^{T} x' J_{\phi_{\theta}}(m_t) \mathrm{Var}(Y_{1,t}) J_{\phi_{\theta}}(m_t)'x \\ &\geq \frac{1}{T} \sum_{t=1}^{T} \gamma_t |J_{\phi_{\theta}}(m_t)'x|^2 \\ &\geq \min_{1 \leq t \leq T} \gamma_t (x' A_T)^2 \\ &\geq \min_{1 \leq t \leq T} \gamma_t \beta_T \\ &\geq \min_{1 \leq t \leq T} \frac{\det \mathrm{Var}(Y_{1,t})}{\mathrm{Tr} \mathrm{Var}(Y_{1,t})} \cdot \frac{\det^2 A_T}{\sum_{1 \leq i,j \leq 2} |A_T(i,j)|^2}. \end{aligned} $$

Let us comment on the previous inequalities. The second and the fourth one are based on the Courant-Fischer theorem. The third one is based on Jensen's inequality. The last inequality is a consequence of the inequality $\lambda \geq \det(B)/\mathrm{Tr}(B)$ for any (non-null) non-negative definite matrix $B$ with an eigenvalue $\lambda$. The result announced follows from the Courant-Fischer theorem and the fact that the assumptions of Theorem 3 guarantee that the coefficients of the matrices $A_T$ and $\mathrm{Var}(Y_{1,t})$ are bounded with respect to $T$ and $t$, respectively.

2. Since $\mu_t = \sigma_t^2$ for the Poisson distribution, using the expressions of the first four moments of the Poisson distribution, we get $\det \mathrm{Var}(Y_{1,t}) = 2\mu_t^3$. The assumption on the function $\mu$ entails that $\inf_{t \geq 1} \mu_t > 0$. Moreover,

$$ A_T = \begin{pmatrix} -2 - \frac{1}{T} \sum_{t=1}^{T} \frac{1}{\mu_t} & \frac{1}{T} \sum_{t=1}^{T} \frac{1}{\mu_t} \\ -\frac{1}{T} \sum_{t=1}^{T} \left(2 + \frac{1}{\mu_t}\right) \log \mu_t & \frac{1}{T} \sum_{t=1}^{T} \frac{\log \mu_t}{\mu_t} \end{pmatrix}. $$

Then

$$ \lim_{T \to \infty} \det(A_T) = 2 \int_0^1 \frac{1}{\mu(u)} du \cdot \int_0^1 \log \mu(u) du - 2 \int_0^1 \frac{\log \mu(u)}{\mu(u)} du. $$

One can show that the previous limit is positive, since from Jensen's inequality applied first to the probability measure $P(du) = 1_{[0,1]}(u)du / (\mu(u) \int_0^1 \mu(v)^{-1}dv)$ and then to the uniform distribution over $[0, 1]$,

$$ \left( \int_0^1 \frac{1}{\mu(u)} du \right)^{-1} \int_0^1 \frac{\log \mu(u)}{\mu(u)} du < - \log \int_0^1 \frac{1}{\mu(u)} du < \int_0^1 \log \mu(u) du. $$

This proves the lemma.

□
**Lemma 7.** Suppose that all the assumptions of Theorem 3 are valid, with or without the additional condition on $T/n$ given in point 1 or point 2. Set $\Sigma_t := \text{Var}(Y_{1,t})$ and

$$ \hat{\Sigma}_t := \frac{1}{n} \sum_{j=1}^{n} Y_{j,t} Y'_{j,t} - \bar{Y}_{n,t} \bar{Y}'_{n,t}, \quad \bar{Y}_{n,t} := \frac{1}{n} \sum_{j=1}^{n} Y_{j,t}, $$

$$ \hat{\Gamma}_T := \frac{1}{T} \sum_{t=1}^{T} J_{\phi_{\hat{\theta}}} (\hat{m}_t) \hat{\Sigma}_t J_{\phi_{\hat{\theta}}} (\hat{m}_t)', $$

$$ \hat{G}_t := (\text{Tr}(\nabla^{(2)}\phi_{1,\hat{\theta}}(\hat{m}_t)\hat{\Sigma}_t), \text{Tr}(\nabla^{(2)}\phi_{2,\hat{\theta}}(\hat{m}_t)\hat{\Sigma}_t)). $$

Then, with $\Gamma_T$ defined in the statement of Lemma 6,

$$ \hat{\Gamma}_T - \Gamma_T = o_p(1). \qquad (26) $$

Moreover, if $C_T := \Gamma_T^{-1/2}$ and $\hat{C}_T := \hat{\Gamma}_T^{-1/2}$, we have

$$ \hat{C}_T - C_T = o_p(1), \qquad (27) $$

$$ \frac{1}{T} \sum_{t=1}^{T} (\hat{G}_t - G_t) = o_p(1). \qquad (28) $$

*Proof.* Under the assumptions of Theorem 3, $\phi_\theta$ and its Jacobian and Hessian matrices are bounded along the sequence $(m_t)_{t \ge 1}$.

We first prove (26). From Lemma 1, it is only necessary to show the convergence on the set

$$ E_n(\epsilon) = \cap_{t=1}^{T} \{ \hat{\mu}_t > \epsilon, \hat{\sigma}_t^2 > \epsilon \}. $$

For such $\epsilon > 0$ and any $\kappa \in (0, 1)$, there exists $d_{\epsilon, \kappa} > 0$ such that for any $z, z' \ge \epsilon$,

$$ |\log z - \log z'| \le d_{\epsilon, \kappa} |z - z'|^{\kappa}, \quad \left| \frac{1}{z} - \frac{1}{z'} \right| \le d_{\epsilon, \kappa} |z - z'|^{\kappa}. $$

Let $0 < \kappa \le 1$ to be fixed later. Using the previous inequalities, one can easily find a constant
$c_{\epsilon, \kappa} > 0$ to satisfy the following inequalities:

$$
\begin{align*}
|\hat{\sigma}_t^2 - \sigma_t^2|^\kappa &\le c_{\epsilon,\kappa} [|\hat{m}_{2,t} - m_{2,t}|^\kappa + |\hat{\mu}_t - \mu_t|^\kappa + |\hat{\mu}_t - \mu_t|^{2\kappa}], \\
|\phi_{1,\hat{\theta}}(\hat{m}_t) - \phi_{1,\theta}(m_t)| &\le c_{\epsilon,\kappa} [|\hat{\sigma}_t^2 - \sigma_t^2|^\kappa + |\hat{\theta}_1 - \theta_1| + |\hat{\theta}_2 - \theta_2| + \{ |\hat{\theta}_2 - \theta_2| + 1 \} \cdot |\hat{\mu}_t - \mu_t|^\kappa], \\
\left|\frac{\partial \phi_{1,\hat{\theta}}(\hat{m}_t) - \partial \phi_{1,\theta}(m_t)}{\partial x}\right| &\le c_{\epsilon,\kappa} [|\hat{\sigma}_t^2 - \sigma_t^2|^\kappa + |\hat{\mu}_t - \mu_t| + |\hat{\theta}_2 - \theta_2| + |\hat{\mu}_t - \mu_t| \cdot |\hat{\sigma}_t^2 - \sigma_t^2|^\kappa], \\
\left|\frac{\partial \phi_{1,\hat{\theta}}(\hat{m}_t) - \partial \phi_{1,\theta}(m_t)}{\partial y}\right| &\le c_{\epsilon,\kappa} |\hat{\sigma}_t^2 - \sigma_t^2|^\kappa, \\
\left|\frac{\partial \phi_{2,\hat{\theta}}(\hat{m}_t) - \partial \phi_{2,\theta}(m_t)}{\partial x}\right| &\le c_{\epsilon,\kappa} [|\hat{\mu}_t - \mu_t|^\kappa + |\phi_{1,\hat{\theta}}(\hat{m}_t) - \phi_{1,\theta}(m_t)| \cdot |\hat{\mu}_t - \mu_t|^\kappa] \\
&\quad + c_{\epsilon,\kappa} \left[ |\phi_{1,\hat{\theta}}(\hat{m}_t) - \phi_{1,\theta}(m_t)| + |\hat{\mu}_t - \mu_t|^\kappa \cdot \left| \frac{\partial \phi_{1,\hat{\theta}}(\hat{m}_t) - \partial \phi_{1,\theta}(m_t)}{\partial x} \right| \right] \\
&\quad + c_{\epsilon,\kappa} \left| \frac{\partial \phi_{1,\hat{\theta}}(\hat{m}_t) - \partial \phi_{1,\theta}(m_t)}{\partial x} \right|, \\
\left|\frac{\partial \phi_{2,\hat{\theta}}(\hat{m}_t) - \partial \phi_{2,\theta}(m_t)}{\partial y}\right| &\le c_{\epsilon,\kappa} [|\hat{\mu}_t - \mu_t|^\kappa + |\hat{\mu}_t - \mu_t|^\kappa \cdot \left| \frac{\partial \phi_{1,\hat{\theta}}(\hat{m}_t) - \partial \phi_{1,\theta}(m_t)}{\partial y} \right|] \\
&\quad + c_{\epsilon,\kappa} \left| \frac{\partial \phi_{1,\hat{\theta}}(\hat{m}_t) - \partial \phi_{1,\theta}(m_t)}{\partial y} \right|.
\end{align*}
$$

Moreover, there exists a constant $c > 0$ such that

$$
\begin{align*}
|\hat{\Sigma}_t(1,2) - \Sigma_t(1,2)| &\le c[|\hat{m}_{3,t} - m_{3,t}| + |\hat{m}_{1,t} - m_{1,t}| \cdot |\hat{m}_{2,t} - m_{2,t}| + |\hat{m}_{1,t} - m_{1,t}| + |\hat{m}_{2,t} - m_{2,t}|], \\
|\hat{\Sigma}_t(2,2) - \Sigma_t(2,2)| &\le c[|\hat{m}_{4,t} - m_{4,t}| + |\hat{m}_{2,t} - m_{2,t}|^2 + |\hat{m}_{2,t} - m_{2,t}|].
\end{align*}
$$

Next, setting $\hat{A}_t = J_{\phi_{\hat{\theta}}}(\hat{m}_t)$ and $A_t = J_{\phi_{\theta}}(m_t)$, we use the decomposition

$$
\begin{align*}
\hat{A}_t \hat{\Sigma} \hat{A}_t - A_t \Sigma_t A_t &= A_t \Sigma_t (\hat{A}_t - A_t) + A_t (\hat{\Sigma}_t - \Sigma_t) A_t + A_t (\hat{\Sigma}_t - \Sigma_t) (\hat{A}_t - A_t) \\
&\quad + (\hat{A}_t - A_t) \Sigma_t A_t + (\hat{A}_t - A_t) \Sigma_t (\hat{A}_t - A_t) + (\hat{A}_t - A_t) (\hat{\Sigma}_t - \Sigma_t) A_t \\
&\quad + (\hat{A}_t - A_t) (\hat{\Sigma}_t - \Sigma_t) (\hat{A}_t - A_t).
\end{align*}
$$

We deduce from the previous bounds that the entries of the difference $\hat{\Gamma}_T - \Gamma_T$ can be decomposed as a sum of terms, each of which is bounded up to a constant possibly depending on $|\hat{\theta}_i - \theta_i|$, $i = 1, 2$, by a partial sum of the form $\frac{1}{T} \sum_{t=1}^{T} \prod_{i=1}^{d} |S_{\alpha_i, t}|^{\kappa_i}$ for some integer $d \ge 1$ and for $i = 1, \dots, 4$, $\kappa_i$ is a positive real number, $\alpha_i \in \{1, \dots, 4\}$ with $\sum_{i=1}^{d} \alpha_i \kappa_i \le 4 + 6\kappa$. Choosing $\kappa \in (0, 1)$ such that $4 + 6\kappa \le k$, one can apply Corollary 3 and use the weak consistency of $\hat{\theta}$ to conclude that $\hat{\Gamma}_T - \Gamma_T = o_P(1)$.

We now prove (27). First we show that

$$
\hat{\Gamma}_T^{-1} - \Gamma_T^{-1} = o_P(1) \tag{29}
$$

using the decomposition

$$
\hat{\Gamma}_T^{-1} - \Gamma_T^{-1} = \hat{\Gamma}_T^{-1} (\Gamma_T - \hat{\Gamma}_T) \Gamma_T^{-1}.
$$
From our assumptions, there exists $K > 0$ such that for $n$ large enough, $|\Gamma_T^{-1}| \le K$. It is then only necessary to show that

$$ \hat{\Gamma}_T^{-1} = O_P(1). \qquad (30) $$

Let $0 < \eta < 1/K$ and set $M := K/(1 - \eta K)$. We have

$$ \mathbb{P}(|\hat{\Gamma}_T^{-1}| > M) \le \mathbb{P}(\sup_{|H|\le\eta} |(\Gamma_T - H)^{-1}| > M) + \mathbb{P}(|\hat{\Gamma}_T - \Gamma_T| > \eta). $$

The second probability on the right side goes to 0 as $n \to \infty$ and the first one is 0 because

$$ (\Gamma_T - H)^{-1} = \sum_{k \ge 0} (\Gamma_T^{-1} H)^k \Gamma_T^{-1} $$

has a norm smaller than $M$. This shows (30) and then (29).

Next, using the Ando-Hemmen inequality (van Hemmen and Ando, 1980, Proposition 3.2), we have for any unitarily invariant norm $|\cdot|,$

$$ \left| \hat{\Gamma}_T^{-1/2} - \Gamma_T^{-1/2} \right| \le \frac{\left| \hat{\Gamma}_T^{-1} - \Gamma_T^{-1} \right|}{\sqrt{\lambda_-(\hat{\Gamma}_T^{-1})} + \sqrt{\lambda_-(\Gamma_T^{-1})}} \le \frac{\left| \hat{\Gamma}_T^{-1} - \Gamma_T^{-1} \right|}{\sqrt{\lambda_-(\Gamma_T^{-1})}}, $$

where $\lambda_-(B)$ denotes the smallest eigenvalue of a symmetric matrix $B$. Of course, the previous inequality is valid for an arbitrary norm $|\cdot|$, up to a constant factor, by equivalence of the norms in a finite dimensional space. If we denote by $\lambda_+(B)$ the largest eigenvalue of a matrix $B$, we have

$$ \underline{\lim} \lambda_-(\Gamma_T^{-1}) = 1/\overline{\lim} \lambda_+(\Gamma_T) > 0 $$

and $\hat{C}_T - C_T = o_P(1)$.

To prove (28), we proceed as for the proof of (26). We first bound the differences

$$ \left| \frac{\partial^2 \phi_{\ell, \hat{\theta}}}{\partial x^i \partial y^j} (\hat{m}_t) - \frac{\partial^2 \phi_{\ell, \theta}}{\partial x^i \partial y^j} (m_t) \right|, \quad \ell = 1, 2, \quad i, j \in \{0, 1, 2\} \text{ such that } i + j = 2. $$

We find that each difference can be bounded by a sum of terms of the form

$$ C_1 |\hat{\theta}_i - \theta_i|, \quad \text{or} \quad C_2 \prod_{k=1}^{d} |S_{\alpha_k, t}|^{\kappa_k} $$

with $C_1, C_2 > 0$, $\alpha_k \in \{1, 2\}$, $\kappa_k \ge 0$ and $\sum_{k=1}^d \alpha_k \kappa_k$ bounded by $2 + 2\kappa$ if $(\ell, i, j) = (1, 2, 0)$, $1 + 2\kappa$ if $(\ell, i, j) = (1, 1, 1)$, $2\kappa$ if $(\ell, i, j) = (1, 0, 2)$, $2 + 3\kappa$ if $(\ell, i, j) = (2, 2, 0)$, $1 + 3\kappa$ for $(\ell, i, j) = (2, 1, 1)$ and $3\kappa$ if $(\ell, i, j) = (2, 0, 2)$. It is then easily seen that the differences $\frac{1}{T} \sum_{t=1}^T (\hat{G}_t - G_t)$ only involve some factors of the previous type with $\sum_{k=1}^d \alpha_k \kappa_k \le 4 + 3\kappa$ where $\kappa > 0$ can be arbitrarily small. We then get the result from Corollary 3, choosing $\kappa$ small enough such that $4 + 3\kappa \le k$. $\square$
# 10 Acknowledgements

This work was funded by CY-AS ("Investissements d'Avenir" ANR-16-IDEX-0008), "EcoDep" PSI-AAP2020-0000000013. P.D. was also funded by the chair FIME: https://fime-lab.org/. J.E.C. thanks Roseanne Benjamin for her help during this work.

# References

Shaul K Bar-Lev. Independent, tough identical results: the class of Tweedie on power variance functions and the class of Bar-Lev and Enis on reproducible natural exponential families. *Int. J. Stat. Probab*, 9(1):30-35, 2020.

Shaul K Bar-Lev and Peter Enis. Reproducibility and natural exponential families with power variance functions. *Annals of Statistics*, 14(4):1507-1522, 1986.

Shaul K. Bar-Lev and Osnat Stramer. Characterizations of natural exponential families with power variance functions by zero regression properties. *Probability theory and related fields*, 76(4):509-522, 1987.

Maurice S. Bartlett. Some notes on insecticide tests in the laboratory and in the field. *Supplement to the Journal of the Royal Statistical Society*, 3(2):185-194, 1936.

Patrick Billingsley. *Convergence of Probability Measures*. John Wiley & Sons, 2013.

C. I. Bliss. Statistical problems in estimating populations of Japanese beetle larvae. *Journal of Economic Entomology*, 34:221-232, 1941. URL https://doi.org/10.1093/jee/34.2.221.

Mark Brown, Joel E. Cohen, and Victor H. de la Peña. Taylor’s law, via ratios, for some distributions with infinite mean. *Journal of Applied Probability*, 54(3):657-669, 2017.

Mark Brown, Joel E. Cohen, Chuan-Fa Tang, and Sheung Chi Phillip Yam. Taylor’s law of fluctuation scaling for semivariances and higher moments of heavy-tailed data. *Proceedings of the National Academy of Sciences*, 118(46):e2108031118, 2021.

S. J. Clark and J. N. Perry. Small sample estimation for Taylor’s power law. *Environmental and Ecological Statistics*, 1:287-302, 1994.

Matthew RD Cobain, Markus Brede, and Clive N Trueman. Taylor’s power law captures the effects of environmental variability on community structure: An example from fishes in the north sea. *Journal of Animal Ecology*, 88(2):290-301, 2019.

Joel E. Cohen. Stochastic population dynamics in a Markovian environment implies Taylor’s power law of fluctuation scaling. *Theoretical Population Biology*, 93:30-37, 2014a. doi: 10.1016/j.tpb.2014.01.001.

Joel E. Cohen. Taylor’s law and abrupt biotic change in a smoothly changing environment. *Theoretical Ecology*, 7(1):77-86, 2014b. doi: 10.1007/s12080-013-0199-z.

Joel E Cohen and Thierry E Huillet. Taylor’s law for some infinitely divisible probability distributions from population models. *Journal of Statistical Physics*, 188(3):33, 2022.
Joel E. Cohen and Takashi Saitoh. Population dynamics, synchrony, and environmental quality of Hokkaido voles lead to temporal and spatial Taylor's laws. *Ecology*, 97(12):3402–3413, 2016.

Joel E. Cohen, Robert Poulin, and Clément Lagrue. Linking parasite populations in hosts to parasite populations in space through Taylor's law and the negative binomial distribution. *Proceedings of the National Academy of Sciences*, 114(1):E47-E56, 2017.

Joel E. Cohen, Richard A. Davis, and Gennady Samorodnitsky. Heavy-tailed distributions, correlations, kurtosis and Taylor's law of fluctuation scaling. *Proceedings of the Royal Society A*, 476 (2244):20200610, 2020.

Joel E. Cohen, Richard A. Davis, and Gennady Samorodnitsky. COVID-19 cases and deaths in the United States follow Taylor's law for heavy-tailed distributions with infinite variance. *Proceedings of the National Academy of Sciences*, 119(38):e2209234119, 2022.

Victor de la Pena, Paul Doukhan, and Yahia Salhi. Taylor's law of fluctuation scaling for semivariances and higher moments of heavy-tailed data. *Journal of Applied Probability*, 59(2):584–607, 2022.

Jérôme Dedecker, Paul Doukhan, Gabriel Lang, León R José Rafael, Sana Louhichi, and Clémentine Prieur. *Weak Dependence*. Springer, 2007.

Paul Doukhan, George Oppenheim, and Murad Taqqu. *Theory and applications of long-range dependence*. Springer Science and Business Media, 2002.

Masami Fujiwara and Joel E. Cohen. Mean and variance of population density and temporal Taylor's law in stochastic stage-structured density-dependent models of exploited fish populations. *Theoretical Ecology*, 8(2):175–186, 2015. doi: 10.1007/s12080-014-0242-8.

Peter Hall and Christopher C. Heyde. *Martingale Limit Theory and its Application*. Academic Press, 2014.

Wassily Hoeffding. Probability inequalities for sums of bounded random variables. *Journal of the American Statistical Association*, 58(301):13–30, 1963.

A. James and M.J. Plank. On fitting power laws to ecological data. *arXiv preprint arXiv:0712.0613*, 2007.

Bent Jørgensen, Clarice GB Demétrio, Erik Kristensen, Gary T Banta, Hans Christian Petersen, and Matthieu Delefosse. Bias-corrected pearson estimating functions for Taylor's power law applied to benthic macrofauna data. *Statistics & probability letters*, 81(7):749–758, 2011.

Joe N. Perry and Ian P. Woiwod. Fitting Taylor's power law. *Oikos*, pages 538–542, 1992.

Gen Sakoda, Hideki Takayasu, and Misako Takayasu. Tracking Poisson parameter for non-stationary discontinuous time series with Taylor's abnormal fluctuation scaling. *Stats*, 2(1):55–69, 2019.

Gennady Samorodnitsky. Long range dependence. *Foundations and Trends in Stochastic Systems*, 1(3):E163–257, 2007.
Peijian Shi, David A. Ratkowsky, Ningtao Wang, Yang Li, Lei Zhao, Gadi VP Reddy, and Bai-Lian Li. Comparison of five methods for parameter estimation under Taylor's power law. *Ecological complexity*, 32:121-130, 2017.

L Roy Taylor, Ian P. Woiwod, and Joe N. Perry. The density-dependence of spatial behaviour and the rarity of randomness. *The Journal of Animal Ecology*, pages 383-406, 1978.

Lionel Roy Taylor. Aggregation, variance and the mean. *Nature*, 189:732-735, 1961.

Robin A. J. Taylor. *Taylor's Power Law: Order and Pattern in Nature*. Academic press, 2019.

M. C. K. Tweedie. An index which distinguishes between some important exponential families. In *Statistics: Applications and New Directions: Proc. Indian Statistical Institute Golden Jubilee International Conference*, pages 579-604, 1984.

Aad W. Van der Vaart. *Asymptotic Statistics*, volume 3. Cambridge university press, 2000.

J. Leo van Hemmen and Tsuneya Ando. An inequality for trace ideals. *Communications in Mathematical Physics*, 76:143-148, 1980.

Meng Xu and Joel E. Cohen. Spatial and temporal autocorrelations affect Taylor's law for US county populations: Descriptive and predictive models. *PLoS One*, 16(1):e0245062, 2021.

Yang Yang, Han Lin Shang, and Joel E. Cohen. Temporal and spatial Taylor's law: Application to Japanese subnational mortality rates. *Journal of the Royal Statistical Society Series A*, 185(4): 1979-2006, 2022.