2. Preliminaries

In this section we present a number of elementary properties of the model under study which will be used extensively in later sections. A secondary goal of the section is to explain why a maximum-likelihood-based estimation approach is, in our context, problematic.

When aiming to estimate the parameters of the on- and off-time distributions, a first possible approach could be based on maximum likelihood. Such an approach requires the evaluation of the probability

\begin{align*}\ell(n_1,\ldots,n_K)\,:\!=\,{\mathbb P}(A_n(1)=n_1,\ldots,A_n(K)=n_K),\end{align*}

with $n_k\in\{0,\ldots,n\}$ for $k=1,\ldots,K$ . We can then maximize (the logarithm of) this $\ell(n_1,\ldots,n_K)$ over the parameter vector, thus yielding a maximum likelihood estimator.

In the specific case where the on- and off-times are both geometrically distributed, the process $\{A_n(k)\}_{k\in{\mathbb N}}$ is Markovian. More specifically, the probability $\ell(n_1,\ldots,n_K)$ allows an explicit evaluation. Indeed, one can write

\begin{align*}\ell(n_1,\ldots,n_K)={\color{black}{\mathbb P}(A_n(1)=n_1)}\prod_{k=2}^{K} {\mathbb P}(A_n(k) = n_k\,|\, A_n(k-1)=n_{k-1}),\end{align*}

where the transition probabilities ${\mathbb P}(A_n(k) = n_k\,|\, A_n(k-1)=n_{k-1})$ allow explicit evaluation. Indeed, conditional on the value of $A_n(k)$ , $A_n(k+1)$ can be written as the sum of two binomial quantities: with $X\sim {\mathbb G}(p)$ and $Y\sim {\mathbb G}(q)$ ,

\begin{align*}A_n(k+1)\sim {\mathbb B}\textrm{in}(A_n(k), 1-p)+{\mathbb B}\textrm{in}(n-A_n(k),q), \end{align*}

with the two quantities on the right-hand side being independent. We conclude that in this case maximum likelihood can be directly applied in order to estimate the parameters p and q.

In all other cases, however, we do not have access to $\ell(n_1,\ldots,n_K)$ , so that maximum likelihood cannot be (directly) applied. It is noted, though, that we can still set up a computational procedure that outputs the joint moment generating function of $(\boldsymbol{1}(1),\ldots,\boldsymbol{1}(K))$ , where $\boldsymbol{1}(k)$ is the indicator function of an arbitrary edge being on at time k; this procedure is presented in Subsection 2.1. As will be pointed out in Remark 2, the joint moment generating function can be used to evaluate the joint distribution of $(\boldsymbol{1}(1),\ldots,\boldsymbol{1}(K))$ (or any subset of these K indicator functions), which is a quantity we will be working with in Subsection 3.3.

2.1. Joint moment generating function

Before moving to the situation in which we start off in stationarity, we first focus on a single edge that just turned on at time 0, and consider the generic random vector $(\boldsymbol{1}(1),\ldots,\boldsymbol{1}(k))$ that records whether the edge is on at times $1,\ldots,K$ . In this context, our objective is to evaluate the underlying moment generating function. In self-evident notation, we wish to calculate

\begin{align*}M_+({\boldsymbol\theta})\,:\!=\,{\mathbb E}_+ \exp\left(\sum_{k=1}^K \theta_k \boldsymbol{1}(k)\right).\end{align*}

Likewise, for an edge that just turned off, our aim is to evaluate, again in self-evident notation,

\begin{align*}M_-({\boldsymbol\theta})\,:\!=\,{\mathbb E}_- \exp\left(\sum_{k=1}^K \theta_k \boldsymbol{1}(k)\right).\end{align*}

In the sequel, the following objects play a pivotal role:

\begin{align*}v_k({\boldsymbol\theta})\,:\!=\,{\mathbb E}_{k,+} \exp\left(\sum_{i=k}^K \theta_i \boldsymbol{1}(i)\right),\:\:\:\:w_k({\boldsymbol\theta})\,:\!=\,{\mathbb E}_{k,-} \exp\left(\sum_{i=k}^K \theta_i \boldsymbol{1}(i)\right),\end{align*}

where the subscript ‘ $k,+$ ’ (‘ $k,-$ ’, respectively) means that the edge just went on (off, respectively) at time k. It is readily checked that these moment generating functions can be evaluated by solving the following system of linear equations:

\begin{align*} v_k({\boldsymbol\theta}) = \sum_{\ell=1}^{K-k} f_{\ell} \,\exp\left(\sum_{i=k}^{k+\ell-1}\theta_i\right)\,w_{k+\ell}({\boldsymbol\theta})+\sum_{\ell=K-k+1}^\infty f_{\ell} \,\exp\left(\sum_{i=k}^{K}\theta_i\right),\end{align*}

and, likewise,

\begin{align*} w_k = \sum_{\ell=1}^{K-k} g_{\ell} \,v_{k+\ell}({\boldsymbol\theta})+\sum_{\ell=K-k+1}^\infty g_{\ell} .\end{align*}

These are 2K linear equations in equally many unknowns. The system is uppertriangular, allowing for efficient recursive solution. Evidently, $M_+({\boldsymbol\theta})= v_1({\boldsymbol\theta})$ and $M_-({\boldsymbol\theta})= w_1({\boldsymbol\theta})$ .

Now we proceed by considering the situation that the on-off process starts in its stationary version. Using the same reasoning as above, we thus have that

\begin{align*} M({\boldsymbol\theta}) = \varrho\,M_+^\textrm{(res)}({\boldsymbol\theta}) +(1-\varrho)\log M_-^\textrm{(res)}({\boldsymbol\theta}),\end{align*}

where

\begin{align*}M_+^\textrm{(res)} ({\boldsymbol\theta})=\sum_{\ell=1}^{K-1} \bar f_{\ell} \,\exp\left(\sum_{i=1}^{\ell}\theta_i\right)\,w_{1+\ell}({\boldsymbol\theta})+\sum_{\ell=K}^\infty \bar f_{\ell} \,\exp\left(\sum_{i=1}^{K}\theta_i\right),\end{align*}

and, likewise,

\begin{align*}M_-^\textrm{(res)} ({\boldsymbol\theta}) = \sum_{\ell=1}^{K-1} \bar g_{\ell} \,v_{1+\ell}({\boldsymbol\theta})+\sum_{\ell=K}^\infty \bar g_{\ell} .\end{align*}

Remark 2. We finish our procedure to generate the joint moment generating function by pointing out that we can, at least in principle, explicitly convert $M({\boldsymbol\theta})$ into the joint distribution of the vector $(\boldsymbol{1}(1),\ldots,\boldsymbol{1}(K))$ . To demonstrate this procedure, consider the case $K = 3$ (which we need in Subsection 3.3). We wish to compute, in self-evident notation, the eight probabilities $p[0,0,0],p[0,0,1],\ldots,p[1,1,1]$ . To this end, we insert all eight elements from $\{\!-\infty,0\}^3$ into $M({\boldsymbol\theta})$ . For instance, $M(\!-\infty,0,0)$ corresponds to the scenario that at the first time epoch the edge is off, while at the other two epochs the edge can be on or off:

\begin{align*}M(\!-\infty,0,0) = p[0,0,0]+p[0,1,0]+p[0,0,1]+p[0,1,1]. \end{align*}

We thus obtain a non-singular system of eight linear equations in the equally many unknowns, which can be solved in the standard fashion. Observe that $M(0,0,0)=1$ corresponds to the requirement that the eight probabilities have to sum to 1. $\Diamond$

2.2. Saddlepoint approximation

In the previous section we concluded that, except in the case of geometric on- and off-times, we cannot explicitly compute the probability $\ell(n_1,\ldots,n_K)$ . One could still set up an approximate maximum likelihood estimation procedure though, following an approach advocated in [Reference Davison, Hautphenne and Kraus14, Reference de Gunst, Hautphenne, Mandjes and Sollie18]. To this end, we first introduce the Legendre transform (or the convex conjugate) of the joint cumulant generating function $\log M({\boldsymbol\theta})$ :

\begin{align*}I(\boldsymbol{n})\,:\!=\,\sup\nolimits_{\boldsymbol\theta}\left({\boldsymbol\theta}^\top\boldsymbol{n}- n \log M({\boldsymbol\theta})\right), \end{align*}

where (in our context) $\boldsymbol{n}\in\{0,\ldots,n\}^K$ . Write ${\boldsymbol\theta}(\boldsymbol{n})$ for the optimizing argument in the definition of $I(\boldsymbol{n})$ for a given vector $\boldsymbol{n}$ (which is unique, as ${\boldsymbol\theta}\mapsto {\boldsymbol\theta}^\top\boldsymbol{n}\!-\log M({\boldsymbol\theta})$ is concave). Then the saddle-point approximation [Reference Butler9, Reference Daniels13], particularly accurate when n is large, states that

\begin{align*}\log \ell(\boldsymbol{n})\equiv \log \ell(n_1,\ldots,n_K) \approx - J(\boldsymbol{n}) \,:\!=\,-\frac{K}{2}\log(2\pi) - \frac{1}{2}\log s(\boldsymbol{n})-\,I(\boldsymbol{n}),\end{align*}

with

\begin{align*}s(\boldsymbol{n}) \,:\!=\,\left. \textrm{det}\left(\frac{\partial^2}{\partial \theta_i\,\partial\theta_j}n\log M({\boldsymbol\theta})\right)\right|_{{\boldsymbol\theta} ={\boldsymbol\theta}(\boldsymbol{n})}.\end{align*}

The idea of ‘approximate maximum likelihood’ is to maximize $-J(\boldsymbol{n})$ over the parameter vector. As a simplification, one often just maximizes the ‘exponential part’ $-I(\boldsymbol{n})$ only. A formal justification of this type of procedure has recently been provided in [Reference Goodman17].

In the context of our estimation problem, the drawback of the approximate maximum likelihood approach lies in the computational effort. For a given value of the parameter vector, evaluation of $J(\boldsymbol{n})$ or $I(\boldsymbol{n})$ requires the maximization of ${\boldsymbol\theta}^\top\boldsymbol{n}\!-\log M({\boldsymbol\theta})$ over the K-dimensional vector ${\boldsymbol\theta}$ , where it is noted that the computational effort required to compute $M({\boldsymbol\theta})$ is (as a consequence of the underlying uppertriangular structure) quadratic in K. It turns out, in particular for larger values of K, that this computation can be rather time-consuming, despite the fact that ${\boldsymbol\theta}\mapsto {\boldsymbol\theta}^\top\boldsymbol{n}\!-\log M({\boldsymbol\theta})$ is concave. On top of this computation, there is the ‘outer loop’: to numerically identify the (approximate) maximum likelihood estimator, these computations have to be done for a sequence of parameter vectors.

Because of the computational challenges that we are facing, we propose to use an estimator based on the method of moments (to be discussed in the next section) which is of intrinsically lower numerical complexity. A second key advantage of the use of the method of moments is that it allows us to work with other subgraph counts than just edges, as explained in Section 4.

3. Method of moments

Recall that our objective is to estimate the parameters of the on- and off-time distributions from observations of the number of edges $A_n(1),\ldots, A_n(K).$ In this section we propose an estimator based on the method of moments. We first define this estimator and then investigate its asymptotic normality.

3.1. Estimator

As a general principle, in the method of moments one needs to have as many ‘moment equations’ as parameters. In our approach the moment equations are based on the quantities, for $\ell\in{\mathbb N}$ ,

\begin{align*}S_{n,0}(k)\,:\!=\,A_n(k)=\sum_{j = 1}^n\boldsymbol{1}_j(k),\:\:\:\:S_{n,\ell}(k)\,:\!=\,A_n(k)A_n(k+\ell)=\left(\sum_{j = 1}^n\boldsymbol{1}_j(k)\right)\left(\sum_{j = 1}^n\boldsymbol{1}_j(k+\ell)\right).\end{align*}

We proceed by providing explicit expressions for the expectations of these quantities, from which the moment equations can be identified. As we have assumed that the process is in equilibrium at time 1, it remains in equilibrium at any time in $1,\ldots,K.$ We thus obtain, with $s_{n,\ell}\,:\!=\,\mathbb{E}\,S_{n,\ell}(k)$ ,

\begin{align*} s_{n,0}&= \sum_{j = 1}^n \mathbb{E}\left[\boldsymbol{1}_j(k)\right] = n\, {\mathbb P}(\boldsymbol{1}_i(k)=1)=n\varrho. \end{align*}

Along the same lines, using the independence between the per-edge processes,

\begin{align*} s_{n,1}&= \mathbb{E} \left[\sum_{j = 1}^n \boldsymbol{1}_j(k) \boldsymbol{1}_j(k+1)\right] + \mathbb{E}\left[\sum_{j = 1}^n\sum_{i \ne j}\boldsymbol{1}_j(k)\boldsymbol{1}_i(k+1) \right]\\ &= n\, {\mathbb P}(\boldsymbol{1}_i(k)=1)\,{\mathbb P}(\boldsymbol{1}_i(k+1)=1\,|\,\boldsymbol{1}_i(k)=1) + (n^2-n)\big( {\mathbb P}(\boldsymbol{1}_i(k)=1)\big)^2\\ &= n \varrho\, (1-\bar f_1) + (n^2-n) \, \varrho^2.\end{align*}

Similar to the computation of $s_{n,1}$ , we find

(1) \begin{align}s_{n,2} &=n \varrho\, \big((1-\bar f_1-\bar f_2) +\bar f_1g_1\big) + (n^2-n) \, \varrho^2,\\s_{n,3}&= n \varrho\, \big((1-\bar f_1-\bar f_2-\bar f_3) +\bar f_1g_2 + \bar f_1 g_1 (1- f_1)+\bar f_2 g_1\big) + (n^2-n) \, \varrho^2.\notag\end{align}

To understand the expression for $s_{n,2}$ , it should be borne in mind that there are two scenarios in which the edge is on at time k and $k+2$ , namely (in self-evident notation) +++ and + $-$ +. Regarding $s_{n,3}$ , a similar reasoning applies but now there are four scenarios: ++++, + $-$ $-$ +, + $-$ ++, ++ $-$ +. Moment equations for $s_{n,\ell}$ with $\ell= 4,5,\ldots$ can be identified in the same fashion (where the enumeration of all possible scenarios becomes increasingly tedious).

Let $\widehat\mu_{n,K}(\ell)$ be an estimator of $s_{n,\ell}$ , where evident choices are (with $\ell\in{\mathbb N}$ )

\begin{align*}\widehat\mu_{n,K}(0)\,:\!=\,\frac{1}{K}\sum_{k=1}^K A_n(k),\:\:\:\widehat\mu_{n,K}(\ell)\,:\!=\,\frac{1}{K-\ell}\sum_{k=1}^{K-\ell} A_n(k)A_n(k+\ell).\end{align*}

In a model with $L\in{\mathbb N}$ parameters, our method of moments estimator of these parameters is a parameter vector solving the L moment equations

\begin{align*}s_{n,\ell}=\widehat\mu_{n,K}(\ell),\:\:\:\:\ell=0,\ldots, L-1. \end{align*}

We proceed by considering a number of special cases, which will return in the numerical experiments presented in Section 5. In general, the ‘on-probability’ $\varrho$ can be estimated by $\widehat{\varrho}_{n,K} = \widehat{\mu}_{n,K}(0)/n$ .

• • In the case that $X\sim{\mathbb G}(p)$ and $Y\sim{\mathbb G}(q)$ , there are two parameters that need to be estimated. The method of moments therefore involves $s_{n,0}$ and $s_{n,1}$ . The associated moment equations reduce to

  (2) \begin{align} s_{n,0}&= n\varrho = n\frac{q}{p+q},\:\:\:\:\: s_{n,1} =n\varrho\,(1-p) +(n^2-n)\varrho^2.\end{align}
  We conclude that
  \begin{align*}\widehat\varrho_{n,K}=\frac{\widehat\mu_{n,K}(0)}{n},\:\:\:\:\:\widehat p_{n,K}=1-\frac{\widehat\mu_{n,K}(1)-(n^2-n)\widehat\varrho_{n,K}^2}{n\widehat\varrho_{n,K}},\end{align*}
  so that our estimators become
  (3) \begin{align}\notag \widehat p_{n,K}&=\frac{\widehat\mu_{n,K}(0)-\widehat\mu_{n,K}(1)+(1-1/n)(\widehat\mu_{n,K}(0))^2}{\widehat\mu_{n,K}(0)},\\\widehat q_{n,K}&=\frac{\widehat\mu_{n,K}(0)-\widehat\mu_{n,K}(1)+(1-1/n)(\widehat\mu_{n,K}(0))^2}{n-\widehat\mu_{n,K}(0)}.\end{align}

• • In the case $X\sim {\mathbb P}\textrm{ar}(1,\alpha)$ and $Y\sim {\mathbb P}\textrm{ar}(1,\beta)$ , two moment equations are needed. Observe that ${\mathbb E}\,X=\zeta(\alpha)$ and ${\mathbb E}\,Y=\zeta(\beta)$ , where $\zeta(\alpha): = \sum_{i = 1}^{\infty}1/i^{\alpha}$ (for $\alpha>1$ ) is the Riemann zeta function. The moment equations thus become

  \begin{align*} s_{n,0} &= n \varrho = n \, \frac{\zeta(\alpha)}{\zeta(\alpha)+\zeta(\beta)} ,\:\:\:\: s_{n,1} = n\varrho\left(1-\frac{1}{\zeta(\alpha)} \right)+ (n^2-n) \varrho^2.\end{align*}
  It follows, after some elementary calculations, that
  \begin{align*}\widehat\varrho_{n,K}=\frac{\widehat\mu_{n,K}(0)}{n},\:\:\:\:\:\widehat \alpha_{n,K}={\zeta^{-1}}\left(\frac{\widehat\mu_{n,K}(0)}{\widehat\mu_{n,K}(0)-\widehat\mu_{n,K}(1)+(1-1/n)(\widehat\mu_{n,K}(0))^2}\right),\end{align*}
  and, similarly,
  \begin{align*} \widehat\beta_{n,K}&={\zeta^{-1}}\left(\frac{n-\widehat\mu_{n,K}(0)}{\widehat\mu_{n,K}(0)-\widehat\mu_{n,K}(1)+(1-1/n)(\widehat\mu_{n,K}(0))^2}\right).\end{align*}

• • In the case that $X\sim{\mathbb W}(1,\alpha)$ and $Y\sim {\mathbb G}\textrm{eo}(q)$ , entailing that there are two parameters, two moment equations are needed. These are

  \begin{align*} s_{n,0} &= n \varrho = n \, \frac{\chi(\alpha)}{\chi(\alpha)+1/q} ,\:\:\:\: s_{n,1} = n\varrho\left(1-\frac{1}{\chi(\alpha)} \right)+ (n^2-n) \varrho^2\,,\end{align*}
  where $\chi(\alpha) \,:\!=\, \sum_{i=1}^{\infty} e^{-(i-1)^\alpha}$ . It follows that $\widehat q_{n,K}$ is given by (3) and
  \begin{align*} & \widehat \alpha_{n,K} = \chi^{-1} \left(\frac{\widehat\mu_{n,K}(0)}{\widehat\mu_{n,K}(0)-\widehat\mu_{n,K}(1)+(1-1/n)(\widehat\mu_{n,K}(0))^2}\right).\end{align*}

• • In the case that $X\sim {\mathbb P}\textrm{ar}(C,\alpha)$ and $Y\sim {\mathbb G}\textrm{eo}(q)$ there are three parameters, so that three moment equations are needed:

  \begin{align*} s_{n,0} &= n \varrho = n \, \frac{\zeta(C,\alpha)}{\zeta(C,\alpha)+1/q} ,\:\:\:\: s_{n,1} = n\varrho\left(1-\frac{1}{\zeta(C,\alpha)} \right)+ (n^2-n) \varrho^2\,, \\ s_{n,2} &= n\varrho\left(1-\frac{1}{\zeta(C,\alpha)}-\frac{(C/(C+1))^{\alpha}}{\zeta(C,\alpha)} + \frac{q}{\zeta(C,\alpha)}\right) + (n^2-n) \varrho^2 \\ &= s_{n,1} + n\varrho\left(\frac{q}{\zeta(C,\alpha)}-\frac{(C/(C+1))^{\alpha}}{\zeta(C,\alpha)} \right) \,,\:\:\:\:\:\:\: \zeta(C,\alpha)\,:\!=\,\sum_{i=1}^\infty \frac{C^\alpha}{(C+(i-1))^\alpha}.\end{align*}
  As it turns out, we again have that $\widehat q_{n,K}$ is given by (3). The values of $\widehat C_{n,K}$ and $\widehat \alpha_{n,K}$ are obtained by solving two equations in two unknowns:
  \begin{align*}\zeta(C,\alpha)&=\frac{\widehat\mu_{n,K}(0)}{\widehat\mu_{n,K}(0)-\widehat\mu_{n,K}(1)+(1-1/n)(\widehat\mu_{n,K}(0))^2}, \\ \left(\frac{C}{C+1}\right)^\alpha &= \widehat q_{n,K}- \frac{\widehat\mu_{n,K}(2)-\widehat\mu_{n,K}(1)}{\widehat\mu_{n,K}(0)-\widehat\mu_{n,K}(1)+(1-1/n)(\widehat\mu_{n,K}(0))^2}.\end{align*}

Remark 3. It is clear that $s_{n,0}$ depends on the distributions of X and Y through ${\mathbb E}\,X$ and ${\mathbb E}\,Y$ only. Interestingly, this holds for $s_{n,1}$ as well: noting that $\bar{f}_1 = 1/{{\mathbb E}\,X}$ , we see that $s_{n,1}$ is also a function of ${\mathbb E}\,X$ and ${\mathbb E}\,Y$ . This means that in instances in which the on- and off-times both have one parameter, one can equivalently estimate ${\mathbb E}\,X$ and ${\mathbb E}\,Y$ and then find the parameter values that are in line with these values.

In this paper we primarily focus on on- and off-times with relatively few parameters. The number of moment equations should clearly match the number of parameters to be estimated. This means that a complication lies in the fact that for larger values of k the covariance between $A_n(1)$ and $A_n(k)$ becomes small, so that $s_{n,k}$ ‘provides little information’ (i.e. $s_{n,k}$ will become increasingly hard to distinguish from $n^2\varrho^2$ as k grows). Consequently, in practical terms, there are evident limits to the number of parameters our method can handle.

In the cases in which both X and Y are characterized through a single parameter, we work with moment equations based on $s_{n,0}$ and $s_{n,1}$ . In principle, however, we could have also set up an estimator based on, say, $s_{n,0}$ and $s_{n,2}$ . We have already pointed out a drawback of that choice: for larger k, $s_{n,k}$ becomes less informative. In addition, as indicated by the involved expression for $s_{n,2}$ given in (1), the moment equations are less convenient to work with. $\Diamond$

3.2. Asymptotic normality for geometric on- and off-times

We proceed by establishing asymptotic normality for the case that $X\sim{\mathbb G}(p)$ and $Y\sim{\mathbb G}(q)$ , applying the delta method; see e.g. [Reference van der Vaart and Wellner28, Section 3.9]. The main idea is that under the proviso that $(\widehat \mu_{n,K}(0), \widehat \mu_{n,K}(1))$ is asymptotically normal (as $K\to\infty$ , that is), then this property is inherited by $(\widehat p_{n,K}, \widehat q_{n,K})$ by virtue of the delta method. The starting point is that we assume the vector $(\widehat \mu_{n,K}(0), \widehat \mu_{n,K}(1))$ satisfies the following convergence in distribution: as $K\to\infty$ ,

(4)

with $\boldsymbol{Z}$ being a bivariate zero-mean normal vector; we justify this claim later in this subsection. At this point one can immediately apply the results from [Reference van der Vaart and Wellner28, Section 3.9], via which one can establish the joint central limit theorem for $(\widehat p_{n,K}, \widehat q_{n,K})$ by evaluating a matrix of partial derivatives, corresponding to the mapping between $(\widehat \mu_{n,K}(0), \widehat \mu_{n,K}(1))$ and $(\widehat p_{n,K}, \widehat q_{n,K})$ . The claim of Proposition 1 later in this section is now a direct consequence of the general result presented in [Reference van der Vaart and Wellner28, Thm. 3.9.4]. However, we choose to also include an intuitive, informal calculation here, as it reveals how properties of $(\widehat \mu_{n,K}(0), \widehat \mu_{n,K}(1))$ carry over to properties of $(\widehat p_{n,K}, \widehat q_{n,K})$ .

The informal reasoning is as follows.

• • The distributional limit (4) can be written equivalently as

  (5)

  
  as $K\to\infty$ . The scaling in (4)–(5), in combination with the definitions of $\widehat \mu_{n,K}(0)$ , $\widehat \mu_{n,K}(1)$ , and $s_{n,0}$ , $s_{n_1}$ suggests that the random vector $\boldsymbol{Z}$ , with mean zero, should have a covariance matrix that is given by the asymptotic rescaled variances and covariances
  (6) \begin{align}v_0 &\,:\!=\,\lim_{K\to\infty}\frac{1}{K}{\mathbb V}\textrm{ar}\left[\sum_{k=1}^K A_n(k)\right],\:\: \notag v_1 \,:\!=\,\lim_{K\to\infty}\frac{1}{K-1}{\mathbb V}\textrm{ar}\left[\sum_{k=1}^{K-1} A_n(k)A_n(k+1)\right],\:\:\\c_{01}&\,:\!=\,\lim_{K\to\infty}\frac{1}{K-1}{\mathbb C}\textrm{ov}\left[\sum_{k=1}^{K} A_n(k), \sum_{k=1}^{K} A_n(k)A_n(k+1)\right].\end{align}
  Proposition 1 confirms this intuition and computes these $v_0$ , $v_1$ , and $c_{01}$ for the case of $X\sim{\mathbb G}(p)$ and $Y\sim{\mathbb G}(q)$ , relying on the definitions in Subsection 3.3 we point out how the three quantities in (6) can be calculated in general, and, in particular, we show that they are finite in a heavy-tailed example.

• • As a second step, based on Taylor approximations of $(\widehat \mu_{n,K}(0), \widehat \mu_{n,K}(1))$ around $(s_{n,0},s_{n,1})$ , we write

  \begin{align*}\widehat p_{n,K}= \frac{s_{n,0}+Z_0/\sqrt{K}-(s_{n,1}+Z_1/\sqrt{K})+(1-1/n)(s_{n,0}+Z_0/\sqrt{K})^2}{s_{n,0}+Z_0/\sqrt{K}}{\color{black}\,+\,o_{\mathbb P}\left(\frac{1}{\sqrt{K}}\right)},\end{align*}
  and an analogous expression for $\widehat q_{n,K}$ . From the moment equations (2) it follows immediately that
  \begin{align*}\frac{s_{n,0}-s_{n,1}+(1-1/n)s_{n,0}^2}{s_{n,0}}=p.\end{align*}
  Upon combining this identity with the familiar linear approximation $(1+x)^{-1} =$ $1 - x +O(x^2)$ as $x\downarrow 0$ , we obtain that
  \begin{align*}\widehat p_{n,K} &= p +\frac{Z_0}{\sqrt{K}}\left(\frac{s_{n,1}}{s_{n,0}^2}+ 1-\frac{1}{n}\right)+ \frac{Z_1}{\sqrt{K}}\left(-\frac{1}{s_{n,0}}\right){\color{black}\,+\,o_{\mathbb P}\left(\frac{1}{\sqrt{K}}\right)}.\end{align*}

• • We thus obtain that $\widehat p_{n,K}$ is asymptotically normally distributed around p, and there is a similar relation for $\widehat q_{n,K}$ : for vectors ${\boldsymbol\gamma}, {\boldsymbol\delta} \in {\mathbb R}^2$ ,

  \begin{align*}\widehat p_{n,K} &= p + \frac{1}{\sqrt{K}}\left(\gamma_{0}Z_0+\gamma_{1}Z_1\right){\color{black}\,+\,o_{\mathbb P}\left(\frac{1}{\sqrt{K}}\right)},\\\widehat q_{n,K} &= q + \frac{1}{\sqrt{K}}\left(\delta_{0}Z_0+\delta_{1}Z_1\right){\color{black}\,+\,o_{\mathbb P}\left(\frac{1}{\sqrt{K}}\right)}.\end{align*}
  It takes a minor computation to verify that $\boldsymbol{\delta} = (q/p)\boldsymbol{\gamma}$ , where the entries of ${\boldsymbol\gamma}$ equal
  (7) \begin{equation} \gamma_0=\frac{s_{n,1}}{s_{n,0}^2}+ 1-\frac{1}{n} =2+\frac{p-q-(p+q)p}{nq},\:\:\:\:\:\gamma_1=-\frac{p+q}{nq}.\end{equation}
  In this way we have identified the parameters pertaining to the limiting distribution of $(\widehat p_{n,K}, \widehat q_{n,K})$ .

The following result presents the asymptotic normality of our estimator in the regime $K\to\infty.$

Proposition 1. Consider the case that $X\sim{\mathbb G}(p)$ and $Y\sim{\mathbb G}(q)$ . Define, with $\gamma_0$ and $\gamma_1$ as given in (7),

(8) \begin{equation} \sigma^2\,:\!=\,\gamma_0^2v_0+2\gamma_0\gamma_1c_{01}+\gamma_1^2v_1,\quad \tau^2\,:\!=\,\sigma^2\,(q/p)^2,\quad \rho\,:\!=\,\sigma^2 (q/p).\end{equation}

Then the following claims apply.

1. (i) The distributional convergence (4) is valid: as $K\to\infty$ ,

   \begin{align*}\sqrt{K}\left(\begin{array}{c}\widehat \mu_{n,K}(0)-s_{n,0}\\\widehat \mu_{n,K}(1)-s_{n,1}\end{array}\right) \overset{\textrm{d}}{\to}{\mathscr N}\left(\left(\begin{array}{c}0\\0\end{array}\right), \left(\begin{array}{cc}v_0&c_{01}\\c_{01}&v_1\end{array}\right)\right),\end{align*}
   with $v_0$ , $v_1$ , and $c_{01}$ as defined in (6).

2. (ii) As $K\to\infty$ ,

   \begin{align*}\sqrt{K}\left(\begin{array}{c}\widehat p_{n,K} -p\\\widehat q_{n,K} -q\end{array}\right)\to {\mathscr N}\left(\left(\begin{array}{c}0\\0\end{array}\right), \left(\begin{array}{cc}\sigma^2&\rho\\\rho&\tau^2\end{array}\right)\right),\end{align*}
   with $\sigma^2$ , $\tau^2$ , and $\rho$ given in (8).

3. (iii) In addition,

   (9) \begin{equation} v_0= n\,\varrho(1-\varrho)\frac{1+\lambda}{1-\lambda}=n\,\frac{pq(2-p-q)}{(p+q)^3}, \end{equation}
   and $v_1$ is given by (A5) and $c_{01}$ by (A7).

Proof. Noting that above we have established part (ii) under the assumption that part (i) applies, we now verify the validity of (4). In the case of $X \sim {\mathbb G}(p)$ , $Y \sim {\mathbb G}(q)$ , we have that $\{A_n(k),A_n(k+1)\}$ is an irreducible Markov chain on the finite state-space $\{0,\ldots,n\}\times\{0,\ldots,n\}$ that we have assumed to be in stationarity. As a consequence, a standard central limit theorem for irreducible finite-state Markov chains in stationarity, states that

\begin{align*}\sqrt{K}\left(\frac{1}{K} \sum_{k=1}^K f(A_n(k), A_n(k+1)) - {\mathbb E} f(A_n(1), A_n(2))\right) \end{align*}

converges in distribution, as $K\to\infty$ , to a centered normal random variable with a variance that depends on the function f chosen, namely

\begin{align*} \sigma^2_f\,:\!=\,{\mathbb V}\textrm{ar} f(A_n(1), A_n(2)) + 2\sum_{k=1}^\infty {\mathbb C}\textrm{ov}\big(f(A_n(1), A_n(2)),f(A_n(k+1), A_n(k+2))\big).\end{align*}

Observe that $v_0$ and $v_1$ from (6) correspond to $\sigma^2_f$ for the cases $f(x,y)=x$ and $f(x,y)=xy$ , respectively. Supposing that we pick $f(x,y)\,:\!=\,\alpha x + \beta xy$ for arbitrary scalars $\alpha,\beta$ , then we find that $\sigma^2_f = \alpha^2\,v_0+\beta^2\, v_1+2\alpha\beta\, c_{01}$ and that

\begin{align*}\sqrt{K} \Big(\alpha\widehat \mu_{n,K}(0)+\beta\widehat \mu_{n,K}(1)-\alpha s_{n,0}-\beta s_{n,1}\Big) \end{align*}

converges to a centered normal random variable, so that with the Cramér–Wold device [Reference Billingsley5, Thm. 29.4] it follows that

\begin{align*}\sqrt{K} \Big(\widehat \mu_{n,K}(0)- s_{n,0},\widehat \mu_{n,K}(1)- s_{n,1}\Big) \end{align*}

converges to a centered bivariate normal random vector. It is noted that the finiteness of $v_0$ , $v_1$ , and $c_{01}$ is guaranteed by part (iii).

We now turn to part (iii), demonstrating how the quantities introduced in (6) can be computed. As we are in the setting of $X\sim{\mathbb G}(p)$ and $Y\sim{\mathbb G}(q)$ , we can exploit the Markovian nature of the process $\{A_n(k)\}_{k\in{\mathbb N}}$ . Define $\lambda\,:\!=\,1-p-q.$ We focus here on the computation of $v_0$ . The quantities $v_1$ and $c_{01}$ can be dealt with in a similar way, but require considerably more intricate calculations; these are provided in Appendix A. We are to evaluate

\begin{align*} {\color{black}v_0} &= {\color{black}{\mathbb V}\mbox{ar}\,A_n(1) + 2 \sum_{k=2}^\infty {\mathbb C}\mbox{ov}\left(A_n(1),A_n(k)\right).} \end{align*}

Now notice that, conditional on $A_n(1)$ , with the two binomial random variables appearing in the right-hand side being independent,

\begin{align*}A_n(k) \sim {\mathbb B}\textrm{in}(A_n(1), r_{k-1}) + {\mathbb B}\textrm{in}(n-A_n(1),s_{k-1}),\end{align*}

where, as follows by an elementary computation,

\begin{align*}r_k&\,:\!=\,{\mathbb P}(\boldsymbol{1}_j(k)=1\,|\,\boldsymbol{1}_j(1)=1)=\varrho +(1-\varrho)\,\lambda^{k-1},\\ s_k&\,:\!=\,{\mathbb P}(\boldsymbol{1}_j(k)=1\,|\,\boldsymbol{1}_j(1)=0)=\varrho -\varrho\,\lambda^{k-1}. \end{align*}

We thus find that

\begin{align*} {\mathbb E}\big[A_n(1)A_n(k)\big] &= {\mathbb E}\big[A_n(1)\big(A_n(1)\,r_{k-1} + (n-A_n(1))\,s_{k-1}\big)\big] = \lambda^{k-1}\,n\varrho(1-\varrho)+n^2\varrho^2, \end{align*}

so that ${\mathbb C}\mbox{ov}\left(A_n(1),A_n(k)\right)=\lambda^{k-1}\,n\varrho(1-\varrho).$ Upon combining the above, and in addition using that ${\mathbb V}\mbox{ar}\,A_n(1)=n\,\varrho(1-\varrho)$ , we obtain that

\begin{align*}v_0 = n\,\varrho(1-\varrho) + 2n\,\varrho(1-\varrho)\sum_{k=2}^{\infty}\,\lambda^{k-1},\end{align*}

which yields the stated expression. As mentioned, the computations for $v_1$ and $c_{01}$ are presented in Appendix A.

Remark 4. In the case $\lambda = 0$ , i.e. $p+q=1$ , it is readily verified that that $\{A_n(k)\}_{k\in{\mathbb N}}$ is a sequence of independent ${\mathbb B}\textrm{in}(n, \varrho)$ random variables. This could suggest that one can only estimate $\varrho$ , and not the individual values of p and q, but this intuition is wrong: the proposed estimators $(\widehat p_{n,K},\widehat q_{n,K})$ are capable of recuperating the true values p and q in an asymptotically normal fashion. The crucial realization is that from the data, as K grows large, it is increasingly picked up that we are dealing with an instance in which $s_{n,1}$ equals $s_{n,0}^2$ . This means that, by (2),

\begin{align*}n^2\varrho^2=n\varrho\,(1-p) +(n^2-n)\varrho^2,\end{align*}

which can be checked to reduce to $p+q=1$ . Being able to estimate $\varrho$ , together with knowing that $p+q=1$ , makes the individual parameters p and q identifiable. Informally speaking, the fact that we learn that there is no correlation between subsequent observations provides helpful information that can be used in the estimation. $\Diamond$

3.3. General recipe for computing (co-)variances

Whereas in the above analysis for the case of geometric on- and off-times we intensively relied on the memoryless property of the geometric distribution, rendering $\{A_n(k)\}_{k\in{\mathbb N}}$ Markov, our next goal is to set up an approach to compute $v_0$ , $v_1$ , and $c_{01}$ for any pair of on- and off-time distributions. In the analysis of the case with geometric on- and off-times (see, in particular, the computations for $v_1$ in Appendix A), we have observed that the most challenging component in the computation is ${\mathbb E}[A_n(1)A_n(2)A_n(k)A_n(k+1)]$ , with the system being in stationarity at time 1. For that reason, we demonstrate how this object can be evaluated for general on- and off-times, leaving it to the reader to adapt this procedure to make it applicable to all objects appearing in $v_0$ , $v_1$ , and $c_{01}$ .

The analysis starts by defining the following probabilities. For $i_2,i_k,i_{k+1}\in\{0,1\}$ ,

\begin{align*}p^-_{i_2,i_k,i_{k+1}} &= {\mathbb P}_\textrm{s}(\boldsymbol{1}(2)=i_2, \boldsymbol{1}(k) = i_k, \boldsymbol{1}(k+1) = i_{k+1}\,|\,\boldsymbol{1}(1) = 0),\\p^+_{i_2,i_k,i_{k+1}} &= {\mathbb P}_\textrm{s}(\boldsymbol{1}(2)=i_2, \boldsymbol{1}(k) = i_k, \boldsymbol{1}(k+1) = i_{k+1}\,|\,\boldsymbol{1}(1) = 1).\end{align*}

Here, the subscript ‘s’ denotes that the process is assumed to be in stationarity at time 1 (which we have assumed throughout). The probabilities $\boldsymbol{p}^\pm\,:\!=\,(p^\pm_{0,0,0},p^\pm_{1,0,0},\ldots, p^\pm_{1,1,1})$ can be determined as in Remark 2.

Now condition on the value of $A_n(1)$ ; say it is $m\in\{0,1,\ldots,n\}$ . Define

\begin{align*}A_n[i_2,i_{k},i_{k+1}] \,:\!=\,\# \big\{\,j\in\{1,\ldots,n\}:\boldsymbol{1}_j(2)=i_2, \boldsymbol{1}_j(k) = i_k, \boldsymbol{1}_j(k+1) = i_{k+1}\big\}. \end{align*}

Then consider the random vector $\boldsymbol{A}_n\equiv (A_n[0,0,0], A_n[0,0,1],\ldots, A_n[1,1,1])\in\{0,\ldots,n\}^8$ . It is evident that $\boldsymbol{A}_n$ , conditional on $A_n(1)=m$ , is the sum of two independent multinomially distributed quantities: in standard notation,

(10) \begin{equation}\boldsymbol{A}_n\sim {\mathbb M}\textrm{ult} \big(m, \boldsymbol{p}^+\big)+{\mathbb M}\textrm{ult} \big(n-m, \boldsymbol{p}^-\big).\end{equation}

Now note that we can write, for a suitably chosen matrix M of dimension $3\times 8$ ,

\begin{align*}\big(\begin{array}{c c c}A_n(2)& A_n(k)& A_n(k+1)\end{array}\big)^\top = M \boldsymbol{A}_n; \end{align*}

for instance, $A_n(2) = A_n[1,0,0]+A_n[1,0,1]+ A_n[1,1,0]+ A_n[1,1,1]$ . Then we write

\begin{align*} {\mathbb E}_\textrm{s}\big[A_n(1)&A_n(2)A_n(k)A_n(k+1)\big]=\sum_{m=0}^n {\mathbb P}(A_n(1)=m)\,m\, {\mathbb E}_\textrm{s}\big[A_n(2)A_n(k)A_n(k+1)\,|\,A_n(1)=m\big]\\ &= \sum_{m=0}^n {\mathbb P}(A_n(1)=m)\,m\,{\mathbb E}_\textrm{s}\big[\big( M \boldsymbol{A}_n\big)_1\big( M \boldsymbol{A}_n\big)_2\big( M \boldsymbol{A}_n\big)_3\,|\,A_n(1)=m\big].\end{align*}

The last step is to plug in (10), and to use that $A_n(1)\sim{\mathbb B}\textrm{in}(n,\varrho).$ Then we end up with a weighted sum of mixed moments (of order up to 4) pertaining to multinomial distributions.

In the light of existing versions of the central limit theorem, it is anticipated that a necessary condition for asymptotic normality is the finiteness of $v_0$ , $v_1$ , and $c_{01}$ . Earlier in this section, this finiteness was verified for the case of geometric on- and off-times by explicitly identifying these three quantities. As demonstrated in Appendix B, one can prove that, in the case where the on-times are ${\mathbb P}\textrm{ar}(C_X,\alpha)$ and the off-times are ${\mathbb P}\textrm{ar}(C_Y,\beta)$ for $\beta>\alpha$ , one needs to have that $\min\{\alpha,\beta\}>2$ , as stated in the following result.

Proposition 2. Consider the case that $X\sim{\mathbb P}\textrm{ar}(C_X,\alpha)$ and $Y\sim{\mathbb P}\textrm{ar}(C_Y,\beta)$ with $C_X, C_Y>0$ . If $\min\{\alpha,\beta\}>2$ , then $v_0$ , $v_1$ , and $c_{01}$ are finite.

Remark 5. Upon inspecting the derivation in Appendix B, it can be seen that the conditions can be relaxed in various respects. In its most general form, one has to require that either ${\mathbb P}(X\geqslant i)$ or ${\mathbb P}(Y\geqslant i)$ has a regularly varying tail [Reference Bingham, Goldie and Teugels6] of index strictly larger than 2, and that the other tail is at most equally heavy. This means, for instance, that $v_0$ , $v_1$ , and $c_{01}$ are finite if $X\sim{\mathbb P}\textrm{ar}(C,\alpha)$ with $\alpha > 2$ and $Y\sim {\mathbb G}(q)$ , or if $X\sim{\mathbb P}\textrm{ar}(C,\alpha)$ with $\alpha > 2$ and $Y\sim {\mathbb W}(\lambda,\beta)$ . $\Diamond$

4. Beyond the number of edges: other subgraph counts

So far, we have set up a parametric inference procedure to estimate the on- and off-time distributions based on observations of the evolution of the total number of edges. However, using essentially the same methodology one can set up a similar procedure if the information available corresponds to other subgraph counts. In this section we demonstrate this for the case that we have access to observations of the number of triangles $T_N(k)$ at times $k\in\{1,\ldots,K\}$ , where $N\in{\mathbb N}$ is the number of vertices of our dynamic Erdős–Rényi graph, and we also present its counterpart for the case that we observe the number of wedges $W_N(k)$ .

Subgraph counts play a crucial role in random graph research. In this context we mention [Reference Chatterjee and Varadhan11], a methodological breakthrough when it comes to the rare-event behavior of static Erdős–Rényi random graphs. Various papers considered the likelihood of certain subgraph counts simultaneously attaining specific values, sometimes giving rise to phase transitions. A few examples of papers in this spirit are [Reference den Hollander, Mandjes, Roccaverde and Starreveld22, Reference Kenyon, Radin, Ren and Sadun23, Reference Radin, Ren and Sadun25]. Regarding the relevance of subgraph counts in various application areas, see e.g. [Reference Gallo, Garlaschelli, Lambiotte, Saracco and Squartini16].

First, we notice that setting up likelihood-based methods would lead to serious complications even in the case of geometrically distributed on- and off-times. To see this, observe that even for geometrically distributed on- and off-times, the process $\{T_N(k)\}_{k\in{\mathbb N}}$ is not Markovian, as one can have rather different configurations that lead to the same number of triangles. As we argue in the remainder of this section, however, one can develop a method of moments.

Denote by $\boldsymbol{1}_{i,j}(k)$ the indicator function that corresponds to the existence of the edge between the vertices i and j at observation $k\in\{1,\ldots,K\}$ , for $i\not= j$ and $i,j\in\{1,\ldots,N\}$ . The first moment equation is

\begin{align*}{\mathbb E}\,T_N(k) = {\mathbb E}\Bigg[\sum_{i_1 < i_2 < i_3}\boldsymbol{1}_{i_1,i_2}(k) \,\boldsymbol{1}_{i_1,i_3}(k)\, \boldsymbol{1}_{i_2,i_3}(k)\Bigg]=\binom{N}{3}\,\varrho^3.\end{align*}

Our next goal is to find a second moment equation. For $k\in\{1,\ldots,K-1\}$ ,

\begin{align*} {\mathbb E}\,T_N(k)&T_N(k+1) \\ &= {\mathbb E}\Bigg[\sum_{i_1 < i_2 < i_3}\boldsymbol{1}_{i_1,i_2}(k) \,\boldsymbol{1}_{i_1,i_3}(k)\, \boldsymbol{1}_{i_2,i_3}(k) \times\! \sum_{j_1 < j_2 < j_3}\boldsymbol{1}_{j_1,j_2}(k+1) \,\boldsymbol{1}_{j_1,j_3}(k+1)\, \boldsymbol{1}_{j_2,j_3}(k+1)\Bigg].\end{align*}

Let $a_m$ , for $m=0,1,2,3$ , be the number of sets $\{\,j_1,j_2,j_3\}$ , for a given set $\{i_1,i_2,i_3\}$ , that have m elements in common. It can be verified that

\begin{align*}a_0=\binom{N-3}{3}=\binom{N}{3}-a_1-a_2-a_3,\:\:\:\:a_1=3\binom{N-3}{2},\:\:\:\: a_2 = 3(N-3),\:\:\:\:a_3 =1.\end{align*}

Observe that (i) if $\{i_1,i_2,i_3\}$ and $\{\,j_1,j_2,j_3\}$ have zero elements or one element in common, then they do not share an edge, (ii) if they have two elements in common, then they share an edge, and (iii) if they are identical, then the corresponding triangles match (in that they share three edges); see Figure 1 for a pictorial illustration. For each instance corresponding to scenario (i), the probability of seeing a triangle at time k as well as $k+1$ is therefore $\varrho^6$ , for each instance corresponding to scenario (ii) it is $\varrho^3\cdot\varrho^2(1-\bar f_1)=\varrho^5(1-\bar f_1)$ , and for each instance corresponding to scenario (iii) it is $\varrho^3(1-\bar f_1)^3.$ It thus follows that

\begin{align*}{\mathbb E}\,T_N(k)T_N(k+1) = \binom{N}{3}\Big(a_0 \varrho^6 + a_1 \varrho^6 + a_2 \varrho^5(1-\bar f_1)+ a_3\varrho^3(1-\bar f_1)^3\Big).\end{align*}

Combining arguments developed earlier in this paper, one can compute ${\mathbb E}\,T_N(k)T_N(k+\ell)$ for any $\ell\in{\mathbb N}$ , so as to obtain additional moment equations (needed if the number of unknown parameters is larger than 2).

Figure 1.

Triangle count, cases $m=1,2$ . The red and black lines represent edges at time k and $k+1$ , respectively. The figure illustrates the elementary fact that (a) if two triangles have one vertex in common ( $m=1$ ), then they do not share an edge, and (b) if they have two vertices in common ( $m=2$ ), then they share one edge.

A similar reasoning can be followed when setting up a method of moments for any type of subgraph count. One could, for instance, consider the number of ‘wedges’ (i.e. groups of three connected vertices)

\begin{align*}W_N(k) = \sum_{i_1\not=i_2<i_3\not=i_1}\boldsymbol{1}_{i_1,i_2}(k) \,\boldsymbol{1}_{i_1,i_3}(k);\end{align*}

here, the index $i_1$ corresponds to the ‘center vertex’, while $i_2$ and $i_3$ are the ‘end edges’. The mean of the wedge count $W_N(k)$ is

\begin{align*}{\mathbb E}\,W_N(k) = {\mathbb E}\Bigg[\sum_{i_1\not=i_2 < i_3\not=i_1}\boldsymbol{1}_{i_1,i_2}(k) \,\boldsymbol{1}_{i_1,i_3}(k)\Bigg]=3 \, \binom{N}{3}\,\varrho^2,\end{align*}

whereas the corresponding second moment equation is

\begin{align*}{\mathbb E}\,W_N(k)W_N(k+1) &= \binom{N}{3}\big(9 a_0 \varrho^4 + 9 a_1 \varrho^4 \nonumber\\&\quad+ {5a_2}\varrho^4 + {4a_2} \varrho^3(1-\bar f_1)+ 3a_3\varrho^2(1-\bar f_1)^2 + 6a_3\varrho^3(1-\bar f_1)\big).\end{align*}

This expression can be understood as follows. To form a wedge, three vertices are needed, and each vertex can be the center node at each of the times, so there are in total nine options.

• • Recall that $a_0$ and $a_1$ relate to the scenarios that $\{i_1,i_2,i_3\}$ and $\{\,j_1,j_2,j_3\}$ have zero elements, respectively one element in common. For the expressions related to $a_0$ or $a_1$ this means that these scenarios do not correspond to an edge being present at times k as well as at $k+1$ .

• • To understand the terms related to $a_2$ , we consider two cases. The first case relates to the configurations in which the two common vertices are both non-center at time k, whereas in the second case one of the vertices is the center vertex at time k. The first case can occur in three ways, where it is directly verified (see the first line of Figure 2(a)) that one cannot have that there is an edge present at time k as well as $k+1$ . The second case can occur in six ways, five of which result in there being no edges in common and one where there is one edge in common (see the middle and bottom rows of Figure 2(a)). Putting the above together, in four out of nine cases there are no edges in common (leading to a contribution $\varrho^4$ ) and in five out of nine cases there is one edge in common (leading to a contribution $\varrho^3(1-\bar f_1)$ ).

• • A similar reasoning applies to the terms related to $a_3$ : if the center vertex at k is the same as at $k+1$ , the edges are the same at both times. As depicted in the first subfigure of Figure 2(b), this gives three options, each of them leading to a contribution $\varrho^2(1-\bar f_1)^2$ ). Otherwise, only one edge is in common. As depicted in the last two subfigures of Figure 2(b), this gives six options, each of them leading to a contribution $\varrho^3(1-\bar f_1)$ ).

Figure 2.

Wedge count, cases $m=2,3$ . (a) $m = 2$ . (Top) If the center node at $k+1$ is not from the vertices that are present at time k, then the two wedges have no edges in common. (Middle) If the center node at $k+1$ is one of the end nodes at k, then the two wedges have one edge in common if the center node at k belongs to the vertices at $k+1$ , and no edge in common otherwise. (Bottom) If the center node at $k+1$ is the center node at k, then there is one edge in common. (b) $m = 3$ . If the center node at $k+1$ is the same as the center node at k, then it is the same wedge; otherwise, there is only one edge in common. The red and black lines represent edges at time k and $k+1$ , respectively.

Figure 2 Long description

The image shows diagrams representing wedge count cases for m equals 2 and 3. In part (a), for m equals 2, it displays three scenarios: the top scenario shows two wedges with no edges in common when the center node at k+1 is not from the vertices present at time k. The middle scenario shows two wedges with one edge in common if the center node at k belongs to the vertices at k+1, and no edge in common otherwise. The bottom scenario shows two wedges with one edge in common if the center node at k+1 is the center node at k. In part (b), for m equals 3, it shows that if the center node at k+1 is the same as the center node at k, it is the same wedge; otherwise, there is only one edge in common. The red and black lines represent edges at time k and k+1, respectively.

5. Numerical experiments

In this section we present numerical experiments by which we empirically assess the performance of our estimation procedure.

In the first series of experiments, we estimate the parameters of the on- and off-times distributions based on the number of edges. We consider the following four instances:

• • $X\sim {\mathbb G}(p)$ with $p=0.3$ , and $Y\sim {\mathbb G}(q)$ with $q=0.8$ .

• • $X\sim {\mathbb P}\textrm{ar}(1,\alpha)$ with $\alpha = 3$ , and $Y\sim {\mathbb P}\textrm{ar}(1,\beta)$ with $\beta = 2.5$ .

• • $X\sim {\mathbb W}(1,\alpha)$ with $\alpha=0.5$ , and $Y\sim {\mathbb G}(q)$ with $q=0.7$ .

• • $X\sim {\mathbb P}\textrm{ar}(C,\alpha)$ with $C=2$ and $\alpha=4$ , and $Y\sim {\mathbb G}(q)$ with $q=0.7$ .

In the first three instances there are two parameters to be estimated, and in the fourth instance three parameters. For each of these four settings the estimators were provided in Subsection 3.1.

We have generated traces of $\{A_n(k)\}$ for $k=1,\ldots,10^5$ and $n=100$ , and estimated all parameters. We performed this procedure $L=1000$ times for each instance; each of the four rows in Figure 3 corresponds to one of the four instances defined above. For example, the two panels in first row of Figure 3, which represent the case where $X\sim {\mathbb G}(p)$ and $Y\sim {\mathbb G}(q)$ , display histograms pertaining to the estimates

\begin{align*}\widehat p^{(1)}_{n,K},\ldots, \widehat p^{(L)}_{n,K},\quad \widehat q^{(1)}_{n,K},\ldots, \widehat q^{(L)}_{n,K}.\end{align*}

In the sequel we denote by

(11) \begin{align}\overline{\widehat p}^{(L)}_{n,K}& \equiv \overline{\widehat p}^{(L)}\,:\!=\, \frac{1}{L}\sum_{\ell=1}^L \widehat p^{(\ell)}_{n,K}, \end{align}

(12) \begin{align} \overline{\widehat q}^{(L)}_{n,K}& \equiv \overline{\widehat q}^{(L)}\,:\!=\, \frac{1}{L}\sum_{\ell=1}^L \widehat q^{(\ell)}_{n,K}\end{align}

the means of these vectors, and by $\widehat\sigma_{p,n,K}^{(L)}\equiv \widehat\sigma_p^{(L)}$ and $\widehat\sigma_{q,n,K}^{(L)}\equiv \widehat\sigma_q^{(L)}$ their standard deviations. The histograms shown in the remaining three rows of Figure 3 should be interpreted in the same manner.

Figure 3.

Estimation from the method of moments when $n=100$ and $K = 10^5$ (results from $L=1000$ experiments). Standard deviations are provided in parentheses.

Figure 3 Long description

The image contains eight histograms representing the distribution of estimated values from a method of moments analysis with 100 experiments and 100,000 iterations. Each histogram shows the frequency distribution of a specific estimated parameter. The x-axis represents the estimated values, while the y-axis represents the frequency of these values. The histograms are color-coded, with some in blue and others in red, indicating different parameters or conditions. Each histogram includes a mean value and a standard deviation in parentheses, providing a summary of the central tendency and variability of the estimates. The distributions vary in shape, with some appearing approximately normal and others skewed. The histograms are labeled with Greek letters and mathematical symbols, indicating the specific parameters being estimated. All values are approximated.

In the experiments reported, the parameters are typically estimated accurately. The exception is the setting of Figure 3(h) where the confidence intervals are wider, potentially a consequence of the fact that in that experiment three (rather than two) parameters are estimated. This means that in that case the moment equations contain $s_{n,0}$ , $s_{n,1}$ , and $s_{n,2}$ (rather than just $s_{n,0}$ and $s_{n,1}$ ), which could have led to the loss of accuracy discussed in Remark Reference Asmussen3. In additional experiments involving Pareto-type distributions (with a fixed number of unknown parameters), not reported in this paper, we have observed that standard deviations increase when the distribution’s shape parameter approaches 2 from above; this behavior is confirmed by the reasoning in Appendix B.

In addition, we generated QQ-plots to assess the asymptotic normality. All of these plots showed nearly straight lines, as desired; see, for example, Figure 4 for those corresponding to the case of geometric on- and off-times. Only in the ‘far tails’ do we observe a significant deviation from the straight line (which will disappear when choosing K larger).

Figure 4.

Q-Q plots of estimates $\overline{\widehat p}^{(L)}$ and $\overline{\widehat q}^{(L)}$ in the ${\mathbb G}(0.3)/{\mathbb G}(0.8)$ case. We picked $n = 100$ and $K =10^5$ (results from $L=1000$ experiments).

Figure 4 Long description

The image contains two Q-Q plots side by side. The left plot shows estimates of p quantiles against standard normal quantiles, with blue plus signs representing data points and a red dashed line indicating the reference line. The right plot shows estimates of q quantiles against standard normal quantiles, similarly represented by blue plus signs and a red dashed reference line. Both plots display a generally linear trend, suggesting that the estimates follow a normal distribution. The data points are closely aligned with the reference line, indicating a good fit. The plots are used to validate the model's assumptions about the distribution of on-times and off-times in a dynamic Erdos-Renyi random graph. The results are based on 1000 experiments with 100 edges and 100,000 iterations.

For the ${\mathbb G}(0.3)/{\mathbb G}(0.8)$ case in Figure 3, we also plotted the probability density functions of a normal distribution in which

• • the respective means are $\overline{\widehat p}^{(L)}$ and $\overline{\widehat q}^{(L)}$ , as given in (11)–(12),

• • the respective variances as given in Proposition 1, with p and q replaced by $\overline{\widehat p}^{(L)}$ and $\overline{\widehat q}^{(L)}$ .

These density functions are in good agreement with the corresponding histograms.

The code to generate the plots is in the checkVariancePQ.m file in the folder Asymptotic in https://github.com/encwang/OnOffEstimation. All the plots in Figures 3 and 5 can be obtained by running the Plot.m file, and the data used in generating the plots in this paper is stored in the .mat file, in each of the folders.

Figure 5.

Estimation from the method of moments based on subgraph counts when $N = 20$ (i.e. $n = 190$ ) and $K = 10^4$ (results from $L=1000$ experiments): (a) estimations based on triangle counts; (b) estimations based on wedge counts. Standard deviations are provided in parentheses.

Figure 5 Long description

The image contains four histograms. The first histogram in the top left corner shows the distribution of estimations based on triangle counts with a mean value of 0.3001 and a standard deviation of 0.0029. The second histogram in the top right corner displays the distribution of estimations based on wedge counts with a mean value of 0.8002 and a standard deviation of 0.0076. The third histogram in the bottom left corner again shows the distribution of estimations based on triangle counts with a mean value of 0.3000 and a standard deviation of 0.0027. The fourth histogram in the bottom right corner shows the distribution of estimations based on wedge counts with a mean value of 0.8001 and a standard deviation of 0.0072. All values are approximated.

In another series of experiments we performed estimations based on triangle and wedge counts. We focused on the case ${\mathbb G}(0.3)/{\mathbb G}(0.8)$ ; we set the number of vertices $N = 20$ , so that $n = N(N-1)/2 = 190$ . In each run we record $\{A_n(k)\}$ for $k = 1,\ldots, 10^4$ , and count the number of triangles (wedges) at each k. The estimation procedure was performed $L=10^3$ times; estimates from these $10^3$ runs have been visualized via histograms in Figure 5. The graphs confirm that the parameters are estimated accurately.