# Statistical Inference in Parametric Preferential Attachment Trees

Fengnan Gao111Fudan University and Shanghai Center for Mathematical Sciences.
Email<EMAIL_ADDRESS>and Aad van der Vaart222TU Delft. Email:
<EMAIL_ADDRESS>
((); (August 16, 2022))

###### Abstract

The preferential attachment (PA) model is a popular way of modelling dynamic
social networks, such as collaboration networks. Assuming that the PA function
takes a parametric form, we propose and study the maximum likelihood estimator
of the parameter. Using a supercritical continuous-time branching process
framework, we prove the almost sure consistency and asymptotic normality of
this estimator. We also provide an estimator that only depends on the final
snapshot of the network and prove its consistency, and its asymptotic
normality under general conditions. We compare the performance of the
estimators to a nonparametric estimator in a small simulation study.

##  1 Introduction and Notation

We study the preferential attachment (PA) model—a dynamic network model which
in our setup evolves from an initial stage consisting of a single node of
degree one (a root node with a loose edge or dead parent) by recursively
adding at each step a single node and edge. The incoming node connects to a
node in the existing network with probability proportional to a non-decreasing
function of its degree. The term preferential attachment reflects that nodes
of higher degrees (“the rich”) inspire more incoming connections (“get
richer”), thus leading to “the-rich-get-richer” effect, or the so-called
Matthew effect.

For a precise description, define $[n]=\\{1,\dots,n\\}$, and denote the nodes
at time $n$ by $\\{v_{i}\\}_{i\in[n]}$. If the corresponding degrees are
$\\{d_{i}(n)\\}_{i\in[n]}$, then the node $v_{n+1}$ connects to the existing
node $v_{i}\in\\{v_{l}\\}_{l\in[n]}$ with probability proportional to
$f(d_{i}(n))$, for a given function $f:\mathbb{N}_{+}\to\mathbb{R}_{+}$, i.e.,
with probability

$\frac{f(d_{i}(n))}{\sum_{j=1}^{n}f(d_{j}(n))}.$

We refer to $f$ as the preferential attachment function, to $f(k)$ as the
preference for a node of degree $k$, and to the denominator in the display as
the total preference at time $n$. After the incoming node $v_{n+1}$ has made
its choice, the scheme repeats itself at time $n+1$ with the set of existing
nodes $\\{v_{i}\\}_{i\in[n+1]}$, their updated degrees
$\\{d_{i}(n+1)\\}_{i\in[n+1]}$, and the incoming node $v_{n+2}$. The model may
evolve to reach any number of nodes, whence we obtain a random graph that
evolves over time.

The PA model received its modern conception and fame in connection to the
prevalence of so-called scale-free networks in the real world, as shown by
Albert Barabási and his co-authors (Barabási and Albert (1999); Barabási,
Albert-László, Albert and Jeong (1999, 2000)), and in many subsequent
scientific studies in different disciplines. A simple version of the PA model
appeared in Barabási and Albert (1999) as a possible explanation of the
emergence of the scale-free property and henceforth the PA model has remained
one of the few _dynamic_ models to produce scale-freeness.

Scale-freeness is usually defined in terms of polynomial decay of the
empirical degree distribution $P_{k}(n)$, which is the proportion of nodes of
degree $k$ at time $n$:

$P_{k}(n)=\frac{1}{n}\sum_{i\in[n]}{1}_{\\{d_{i}(n)=k\\}}=:\frac{1}{n}N_{k}(n).$

In the case that the PA function $f$ is affine with $f(k)=k+{\alpha}$ with
$\alpha>-1$, it is known that $P_{k}(n)\rightarrow p_{k}$ almost surely, as
$n\rightarrow\infty$, for any fixed $k$, where the limit
$(p_{k})_{k=1}^{\infty}$ is a proper probability distribution with atoms
proportional to $k^{-(3+{\alpha})}$, as $k\rightarrow\infty$ (up to a slow-
varying factor, see Móri (2002); van der Hofstad (2017)). The Barabási–Albert
model is the special case with $f(k)=k$ and gives decrease proportional to
$k^{-3}$.

Two different scenarios arise when we no longer restrict attention to the
affine PA functions—superlinear and sublinear. Roughly speaking, in the
superlinear case the PA function $f$ grows faster than any linear function and
the resulting PA tree looks like a star with one dominating node (with high
probability, the first or the oldest one) connecting to virtually almost every
other node, see Oliveira and Spencer (2005) and references therein. In the
sublinear case $f$ grows more slowly than any linear function, yielding more
interesting PA trees. The slower growth of strictly sublinear $f$ yields less
preference towards high-degree nodes, rendering the emergence of high-degree
nodes less likely, and leading to limiting degree distributions with in
general lighter tails than power laws, which have fine and subtle details.

In this paper we consider the statistical estimation of the PA function from
an observed network. We adopt a parametric specification of a sublinear PA
function and consider statistical inference on the parameter. A prototype of
such a model is $f(k)=(k+\alpha)^{\beta}$, for parameters ${\alpha}$ and
$\beta\leq 1$. This includes the submodel $f(k)=k^{\beta}$, for $\beta\leq 1$,
considered in (Barabási, 2016, Section 5.8), and the affine linear model
$f(k)=k+\alpha$. The latter model was considered in Gao and van der Vaart
(2017) by relatively direct arguments and martingale methods, which break down
for more general models. The main contribution of the present paper is to
exploit the framework of supercritical Malthusian branching processes due to
Jagers (1975) and Nerman (1981), and applied to derive the limiting degree
distribution in PA models by Rudas, Tóth and Valkó (2007), to analyze general
parametric models. The branching process framework allows to study the
likelihood function and to prove the consistency and asymptotic normality of
the maximum likelihood estimator.

We also propose a Wald-type test to test the null hypothesis that the PA
function is affine, thus allowing a test for scale-freeness versus a well-
defined alternative. In the case that only a final snapshot and not the
evolution history is observed, we propose a history-free remedy to the
likelihood function and obtain a pseudo maximum likelihood estimator, which is
shown to be consistent in general and asymptotically normal under general
conditions, and those conditions are verified for the case where the PA
function becomes constant after reaching a certain degree.

An alternative to parametric estimation is the nonparametric approach of Gao
et al. (2017), who introduce an empirical estimator and show that this is
consistent for general PA functions $f$. It is unknown whether this empirical
estimator is asymptotically normal. We show by simulation that the maximum
likelihood estimator of the present paper is significantly more efficient if
the parametric model is correctly specified.

###  1.1 Outline of the paper

The paper is organized as follows. We present the likelihood and the maximum
likelihood estimator in Section 2. Section 3 introduces supercritical
Malthusian branching processes and deduce the results needed for the present
paper (with a more formal introduction in the appendix, Section A). In Section
4, we prove that the maximum likelihood estimator is consistent and in Section
5 that it is asymptotically normal. To overcome the problem of relying on the
entire evolution history of the network, we propose in Section 6 a pseudo
maximum likelihood estimator, which depends only on the final snapshot of the
tree, and give conditions for its asymptotic normality. Section 7 gives a new
perspective on the empirical estimator, presenting this as a nonparametric
version of the pseudo maximum likelihood estimator. In Section 8 we present
simulation results that demonstrate the performance of the maximum likelihood
estimator and pseudo maximum likelihood estimator, and compare this to the
empirical estimator. Section 9 collects all proofs to the results. In a second
section B of the appendix we prove the asymptotic normality of the empirical
degree distribution in a new case of special interest, and give a new proof in
the affine case.

###  1.2 Notation

Write $\mathbb{N}_{+}$ for the set of positive natural numbers
$\\{1,2,\dots,\\}$ and $\mathbb{N}$ for the set of natural numbers including
zero. Write $\mathbb{R}_{+}=(0,\infty)$. For a sequence
$(a_{k})_{k=1}^{\infty}$, define $a_{>k}=\sum_{l>k}a_{l}$. Let $N_{k}(t)$ be
the number of nodes of degree $k$ in the network at time $t$, and
$P_{k}(t)=N_{k}(t)/t$ the proportion of such nodes. For a given function
$h:\mathbb{N}_{+}\to\mathbb{R}$ set
$S_{h}(t)=\sum_{k=1}^{\infty}h(k)N_{k}(t)$. The superscript (0) stresses that
a quantity, such as the limiting degree distribution
$p_{k}^{\scriptscriptstyle(0)}$, is considered under the true parameter
$\theta_{0}$. For a vector $v=(v_{1},\dots,v_{d})^{T}\in\mathbb{R}^{d}$, let
$\|v\|:=\|v\|_{\infty}=\max_{1\leq i\leq d}|v_{i}|$; and for a matrix
$A\in\mathbb{R}^{m\times n}$, set $\|A\|:=\max_{i,j}|A_{ij}|$. Define the
diagonal matrix $\mathop{\rm diag\,}\nolimits(a_{1},\dots,a_{d})$ by
$\bigl{(}\mathop{\rm
diag\,}\nolimits(a_{1},\dots,a_{d})\bigr{)}_{ij}={1}_{\\{i=j\\}}a_{i}$. Define
$u_{-}=\max(0,-u)$ for $u\in\mathbb{R}$. We write $a\wedge b=\min(a,b)$ for
$a,b\in\mathbb{R}$. When $cb_{n}\leq a_{n}\leq Cb_{n}$ for some positive
constants $c$ and $C$, we write $a_{n}\asymp b_{n}$.

###  1.3 Model

Throughout the paper $\\{f_{\theta}:\theta\in\Theta\\}$ is a collection of
non-decreasing functions $f_{\theta}:\mathbb{N}_{+}\to\mathbb{R}_{+}$ indexed
by a subset $\Theta\subset\mathbb{R}^{d}$. It is assumed that every element of
the family satisfies one of the two possibilities:

  1. (i)

$f_{\theta}(k)\leq Ck^{\beta}$ for every $k\in\mathbb{N}$ for some positive
constants $C$ and $\beta<1$.

  2. (ii)

$f_{\theta}(k)=k+\alpha$ for every $k\in\mathbb{N}$, for some ${\alpha}>-1$.

It is also assumed that ${\theta}\mapsto f_{\theta}(k)$ is twice continuously
differentiable with derivatives denoted by $\dot{f}_{\theta}(k)$ and
$\ddot{f}_{\theta}(k)$, which are a vector (gradient) and a matrix (Hessian)
in the case of a multidimensional parameter. Differentiation with respect to
the parameter $\theta$ is also denoted by a dot in general.

##  2 Construction of the Maximum Likelihood Estimator

Let $N_{k}(t)$ be the number of nodes of degree $k\in\mathbb{N}$ in the graph
with nodes $v_{1},\ldots,v_{t}$. Initially there is a single node $v_{1}$ with
degree one and hence $N_{1}(1)=1$, and we can set $N_{k}(1)=0$, for every
$k\geq 2$, to define a complete degree sequence. If $D_{t}$ denotes the degree
of the node to which the node $v_{t}$ is attached, then, for $t\geq 2$,

$\displaystyle
N_{k}(t)=N_{k}(t-1)+{1}_{\\{D_{t}=k-1\\}}-{1}_{\\{D_{t}=k\\}}+{1}_{\\{k=1\\}}.$
(2.1)

The random graph evolves as a Markov process, and hence the likelihood
factorizes as the conditional likelihoods of the new node given the current
tree. Any of the existing nodes may be chosen to attach the new node, but only
the degree $D_{t}$ of this node is important for the value of the likelihood,
which takes the form

$L_{n}(f_{\theta})=\prod_{t=2}^{n}\frac{f_{\theta}(D_{t})N_{D_{t}}(t-1)}{S_{f_{\theta}}(t-1)},$

where the norming “constant”
$S_{f_{\theta}}(t-1)=\sum_{k=1}^{\infty}f_{\theta}(k)N_{k}(t-1)$ is the total
preference in the graph with nodes $v_{1},\ldots,v_{t-1}$ given the PA
function $f_{\theta}$. The total preference can be computed recursively by the
rule $S_{f}(t)=S_{f}(t-1)+f(D_{t}+1)-f(D_{t})+f(1)$, for $t\geq 2$, with the
initialization $S_{f}(1)=f(1)$. In particular, the total preference at stage
$t$ can be expressed in the degree sequence up to time $t-1$, and the full
likelihood depends on the data only through $D^{(n)}=(D_{t})_{t=2}^{n}$. The
normalized log-likelihood up to the term $\sum_{t=2}^{n}\log N_{D_{t}}(t-1)$
is given by

$\displaystyle\iota_{n}(f_{\theta})$
$\displaystyle=\frac{1}{n}\sum_{k=1}^{\infty}\log
f_{\theta}(k)\sum_{t=2}^{n}{1}_{\\{D_{t}=k\\}}-\frac{1}{n}\sum_{t=2}^{n}\log
S_{f_{\theta}}(t-1)$ (2.2) $\displaystyle=\sum_{k=1}^{\infty}\log
f_{\theta}(k)P_{>k}(n)-\frac{1}{n}\sum_{t=2}^{n}\log S_{f_{\theta}}(t-1).$

In the last step we use the identity
$\sum_{t=2}^{n}{1}_{\\{D_{t}=k\\}}=N_{>k}(n)$, which is essentially (Gao et
al., 2017, Lemma 1) and results from the fact that any node of degree strictly
large $k$ in the tree at time $n$ must have been chosen for attachment while
it had degree $k$, exactly once up until this time (namely when its degree
went up from $k$ to $k+1$).

The derivative of the $\log$-likelihood is

$\displaystyle\dot{\ell}_{n}(f_{\theta})=\frac{\partial}{\partial\theta}\log
L_{n}(f_{\theta})$
$\displaystyle=\sum_{t=2}^{n}\left[\frac{\dot{f}_{\theta}}{f_{\theta}}(D_{t})-\frac{S_{\dot{f}_{\theta}}(t-1)}{S_{f_{\theta}}(t-1)}\right]$
(2.3)
$\displaystyle=\sum_{t=2}^{n}\left[\frac{\dot{f}_{\theta}}{f_{\theta}}(D_{t})-\mathbb{E}_{{\theta}}\Bigl{[}\frac{\dot{f}_{\theta}}{f_{\theta}}(D_{t})\Big{|}\mathcal{F}_{t-1}\Bigr{]}\right],$

where $(\mathcal{F}_{t})_{t\geq 1}$ is the filtration generated by the
stochastic process of the graph’s evolution. The last expression follows from
the fact that
$\mathbb{P}(D_{t}=k|\mathcal{F}_{t-1})=f_{\theta}(k)N_{k}(t-1)/S_{f_{\theta}}(t-1)$,
and shows that the score process is a martingale under the true parameter
$\theta$, as usual, which can also be seen by readily verifying
$\mathbb{E}[\dot{\ell}_{n}(f_{\theta})\mid\mathcal{F}_{n-1}]=\dot{\ell}_{n-1}(f_{\theta}).$
Dividing this martingale by the number $n$ of nodes in the network and
rewriting as in (2.2) gives

$\displaystyle\dot{\iota}_{n}(f_{\theta})=\frac{1}{n}\dot{\ell}_{n}(f_{\theta})$
$\displaystyle=\sum_{k=1}^{\infty}\frac{\dot{f}_{\theta}}{f_{\theta}}(k)P_{>k}(n)-\frac{1}{n}\sum_{t=2}^{n}\frac{S_{\dot{f}_{\theta}}(t-1)}{S_{f_{\theta}}(t-1)}$
(2.4)
$\displaystyle=\sum_{k=1}^{\infty}\frac{\dot{f}_{\theta}}{f_{\theta}}(k)P_{>k}(n)-\frac{1}{n}\sum_{t=2}^{n}\frac{\sum_{k=1}^{\infty}P_{k}(t-1)\dot{f}_{\theta}(k)}{\sum_{k=1}^{\infty}P_{k}(t-1)f_{\theta}(k)}.$

The maximum likelihood estimator $\hat{\theta}_{n}$ can be defined to be
either the maximizer of the log-likelihood (2.3) or a solution to the equation
$\dot{\iota}_{n}(f_{\theta})=0$.

To understand the behavior of the maximum likelihood estimator, we need to
study the rescaled log-likelihood in (2.2) or its derivative (2.4). However,
the quantities in these equations are anything but easy—they are functionals
of the entire evolution history of a complex Markov process, where the
influence of the past persists in the likelihood. Simple and straightforward
approaches such as martingale methods (cf. (van der Hofstad, 2017, Chapter 8))
are no longer meaningful. Instead, we employ the theory of supercritical
Malthusian branching processes, which will be introduced in the section 3.

Adapting such a powerful framework, we will be able to assert the utility of
the maximum likelihood estimator by showing its consistency and asymptotic
normality, respectively. In particular, in both Theorems 4.1 (from the
viewpoint of M-estimator) and 4.6 (from the viewpoint of Z-estimator), we show
under some mild assumptions on the parametric family that
$\hat{\theta}\to\theta_{0}$ in an almost sure sense as the number of nodes
$n\to\infty$. Theorem 5.1 illuminates that
$\sqrt{n}(\hat{\theta}-\theta_{0})\rightsquigarrow
N(0,\sigma_{f_{\theta_{0}}}^{2})$ for some properly defined
$\sigma_{f_{\theta_{0}}}^{2}$ such that constructing confidence sets and
testing particular well-specified hypotheses are possible.

##  3 The continuous random tree model

Before stating the main results formally, it is necessary to adopt a
continuous-time framework, where nodes are added after exponentially
distributed waiting times. To set this up we equip every node $v$ with a pure
birth process $\xi_{v}$ whose events correspond to new, future nodes being
attached to this particular node and which in calendar time starts at its own
birth, i.e. when it is added to the tree. These birth processes are i.i.d.
across nodes and a typical birth process $\xi=(\xi(t):t\geq 0)$ has birth rate
equal to $f\bigl{(}\xi(t)+1\bigr{)}$. Thus, $\xi$ is a continuous-time Markov
process with state space $\mathbb{N}_{+}$, initial value $\xi(0)=0$, and with
the only possible transitions stepwise increases $k-1\rightarrow k$,
determined by

$\mathbb{P}\bigl{(}\xi(t+dt)=k\mid\xi(t)=k-1\bigr{)}=f(k)\,dt+o(dt).$ (3.1)

Every birth corresponds to a new node attached to the existing tree at the
node whose birth process produced the event. A node whose birth process has
had $k-1$ births will have $k-1$ children and one parent in the tree and hence
possess degree $k$. It will produce a new child with rate $f(k)$, explaining
the right side of the display. At a given calendar time $t$ every node in the
current tree $\Upsilon_{t}$ will have a corresponding, active birth process.
The total rate of all active birth processes will be
$S(t)=\sum_{v\in\Upsilon_{t}}f(d_{v}(t))$, where $d_{v}(t)$ is the degree of
$v\in\Upsilon_{t}$, and a new node will be attached to node $v\in\Upsilon_{t}$
with probability $f(d_{v}(t))/S(t)$ after an exponential waiting time with
mean $1/S(t)$. The process $\Upsilon_{t}$ starts at time $0$ with a single
node that is understood to have degree 1 and hence the first birth will be
after an exponential time with mean $f(1)$. For a more formal setup in the
language of general branching processes, see Section A.

Thus, we obtain a continuous-time branching process $\Upsilon_{t}$ that
contains the discrete-time process as a skeleton. We define $T_{t}$ as the
total number of births in the continuous process up until time $t$, and
$\tau_{n}=\inf\\{t>0:T_{t}\geq n-1\\}$, for $n=1,2,\dots$, as the time of the
$n$th birth (where $\tau_{1}=0$). When evaluated at the stopping times
${\tau}_{1},{\tau}_{2},\ldots$, the continuous-time process gives a sequence
of trees $\Upsilon_{\tau_{1}},\Upsilon_{\tau_{2}},\ldots$, that is equivalent
to the PA model.

The advantage of the continuous-time setup is that the results on branching
processes become more straightforward. To every node we may attach besides a
birth process a second continuous-time process, called a characteristic, also
starting at the birth of the node. Just as the birth processes, the
characteristics are assumed identically distributed, and the birth process and
characteristic attached to a single node may be dependent, but every time a
new node is added, a pair of a birth process and a characteristic are created
that evolve independently of the processes attached to the nodes that appeared
earlier in the history of the tree. For a given characteristic $\varphi$ we
consider the process $(Z_{t}^{\varphi}:t\geq 0)$ given by

$Z_{t}^{\varphi}=\sum_{v\in\Upsilon_{t}}\varphi_{v}(t-\sigma_{v}).$

Here $\sigma_{v}$ is the calendar time at which node $v$ is added to the tree,
so that $t-\sigma_{v}$ is the lifetime of the node since its birth. The
characteristic of node $v$, denoted by $\varphi_{v}$, has $\varphi_{v}(t)$
interpreted as its value at age $t$ and hence $\varphi_{v}(t-\sigma_{v})$ as
its value at calendar time $t$. The variable $Z_{t}^{\varphi}$ gives the sum
of the characteristics of all individuals in the tree at time $t$. In the
supercritical case the processes $Z_{t}^{\varphi}$ grow exponentially in time
at a rate $\mathrm{e}^{{\lambda^{*}}t}$, where ${\lambda^{*}}$ is the so-
called Malthusian parameter, and $\mathrm{e}^{-{\lambda^{*}}t}Z_{t}^{\varphi}$
tends to a (random) limit as $t\rightarrow\infty$. We shall employ these limit
theorems with appropriate choices of characteristics to derive the asymptotics
of the likelihood function.

A key element is the Laplace transform of the reproduction function
$\mu(t)=\mathbb{E}[\xi(t)]$, the mean number of births of a single node at age
$t$, which in our case can be expressed in the PA function as (see Rudas, Tóth
and Valkó (2007) or the proof in Section 9)

$\rho_{f}(\lambda):=\int_{0}^{\infty}\mathrm{e}^{-\lambda
t}\,\mu(dt)=\sum_{l=1}^{\infty}\prod_{k=1}^{l}\frac{f(k)}{\lambda+f(k)}.$
(3.2)

The function $\rho_{f}$ is convex and decreasing on its domain (the set where
it is finite), which is an interval $(\underline{\lambda},\infty)$ or
$[\underline{\lambda},\infty)$ in the positive half line, and tends to zero as
$\lambda\rightarrow\infty$. The Malthusian parameter ${\lambda^{*}}$ is the
solution of the equation $\rho_{f}({\lambda^{*}})=1$.

In the case of a strictly sublinear PA function $f$, we have
$\underline{{\lambda}}=0$ and $\rho_{f}(\lambda)\uparrow\infty$ as
$\lambda\downarrow 0$, while for the PA function $f(k)=k+\alpha$ we have
$\underline{\lambda}=1$ and the exact form
$\rho_{f}(\lambda)=(1+\alpha)/(\lambda-1)$ is known (see Gao et al. (2017);
Rudas, Tóth and Valkó (2007) or the proof of Proposition 3.1). In both cases
the range of $\rho_{f}$ contains the point $1$ as an interior point and the
Malthusian parameter exists. If the sublinearity assumption is violated in the
sense that neither of the two sublinear conditions holds, the existence of the
Malthusian parameter is not guaranteed, and the associated branching process
could explode or behave irregularly, whence the branching process framework
and the ensuing may results fail.

Furthermore, under the sublinear assumptions of the PA function, we have the
following limit theorem.

######  Proposition 3.1.

Suppose that the range of $\rho_{f}$ contains an open neighborhood of 1 and
Let $\varphi_{1}$ and $\varphi_{2}$ be monotone increasing characteristics
such that, for some constants $C>0$ and $\gamma\geq 0$ and every $t>0$,

$\varphi_{i}(t)\leq C\,\xi(t)^{2},\qquad i=1,2,\quad\text{almost surely}.$
(3.3)

If $f:\mathbb{N}_{+}\to\mathbb{R}_{+}$ is monotone with $f(k)\leq Ck^{\beta}$
for some constants $C$ and ${\beta}<1$, then, as $t\rightarrow\infty$,

$\frac{Z_{t}^{\varphi_{1}}}{Z_{t}^{\varphi_{2}}}\xrightarrow{\text{a.s.}}\frac{\int_{0}^{\infty}\mathrm{e}^{-{\lambda^{*}}t}\mathbb{E}[\varphi_{1}(t)]\,dt}{\int_{0}^{\infty}\mathrm{e}^{-{\lambda^{*}}t}\mathbb{E}[\varphi_{2}(t)]\,dt}.$
(3.4)

The same is true if $f(k)=k+\alpha$, for some ${\alpha}>-1$, provided that for
some $r<2+\alpha$ with $r\leq 2$,

$\varphi_{i}(t)\leq C\,\xi(t)^{r},\qquad i=1,2,\quad\text{almost surely}.$
(3.5)

For a given PA function $f$ with Malthusian parameter ${\lambda^{*}}$, define

$p_{k}=\frac{{\lambda^{*}}}{{\lambda^{*}}+f(k)}\prod_{j=1}^{k-1}\frac{f(j)}{{\lambda^{*}}+f(j)},\qquad
k\in\mathbb{N}_{+},$ (3.6)

where the empty product is defined to be $1$, so that
$p_{1}={\lambda^{*}}/\bigl{(}{\lambda^{*}}+f(1)\bigr{)}$. The definition of
${\lambda^{*}}$ may be used to show that $(p_{k})_{k=1}^{\infty}$ is a
probability distribution on $\mathbb{N}_{+}$. In fact, it is the limit of the
empirical degree distribution $(P_{k}(n))_{k=1}^{\infty}$ of the PA network,
as shown by Rudas, Tóth and Valkó (2007). More generally, we have the
following limiting result.

######  Corollary 3.2.

If $f:\mathbb{N}_{+}\to\mathbb{R}_{+}$ is monotone increasing and satisfies
$f(k)\leq Ck^{\beta}$, for some $\beta<1$, and
$h:\mathbb{N}_{+}\rightarrow\mathbb{R}_{+}$ satisfies $h(k)\leq Ck^{2}$, for a
constant $C$, and every $k$, then the empirical degrees $P_{k}(n)$ in the
model with PA function $f$ satisfy, as $n\rightarrow\infty$,

$\sum_{k=1}^{\infty}h(k)P_{k}(n)\xrightarrow{\text{a.s.}}\sum_{k=1}^{\infty}h(k)p_{k}.$
(3.7)

The same is true for the PA function given by $f(k)=k+{\alpha}$, for some
${\alpha}>-1$, and every $h:\mathbb{N}_{+}\rightarrow\mathbb{R}_{+}$
satisfying $h(k)\leq Ck^{r}$, for some $r<2+\alpha$ with $r\leq 2$.

Choosing $h$ equal to the indicator of the set $\\{k\\}$ for a given $k$, we
recover the convergence $P_{k}(n)\rightarrow p_{k}$ of the empirical degrees
to the limit $p_{k}$, first obtained in Rudas, Tóth and Valkó (2007).

It is worth noting that the tail of the distribution $(p_{k})_{k=1}^{\infty}$
(of which the dependence on the PA function $f$ is suppressed from the
notation), as $k\rightarrow\infty$ is heaviest among sublinear $f$ when $f$ is
affine with $f:k\mapsto k+\alpha$, and it corresponds to the limiting
(asymptotic) power law with exponent $3+\alpha$.

The following lemma records two useful identities, which readily follow from
the definition (3.6) of $p_{k}$, and the definition of ${\lambda^{*}}$ in
terms of the Laplace transform (3.2).

######  Lemma 3.3.

Suppose $(p_{k})_{k=1}^{\infty}$ is the limiting degree distribution specified
in (3.6) for the PA function $f$ with Malthusian parameter ${\lambda^{*}}$.
Then ${\lambda^{*}}=\sum_{j=1}^{\infty}f(j)p_{j}$ and, for all $k\geq 1$,

$\displaystyle p_{>k}$
$\displaystyle=\frac{f(k)p_{k}}{\sum_{j=1}^{\infty}f(j)p_{j}}.$ (3.8)

##  4 Consistency

The following theorem shows that the maximum likelihood estimator in the model
introduced in Section 1.3 is consistent. We assume that the true parameter
${\theta}_{0}\in\Theta$ is identifiable in the model
$\\{f_{\theta}:\theta\in\Theta\\}$ in the sense that
$f_{\theta}(k)=cf_{{\theta}_{0}}(k)$ for every $k\in\mathbb{N}_{+}$ and some
constant $c$ if and only if ${\theta}={\theta}_{0}$.

######  Theorem 4.1.

In the model $\\{f_{\theta}:\theta\in\Theta\\}$ stated in Section 1.3 with
compact parameter space $\Theta\subset\mathbb{R}^{d}$ and PA functions
satisfying $f_{\theta}(k)\leq Ck^{\beta}$ for some constants $C$ and $\beta<1$
or $f_{\theta}(k)=k+{\alpha}$ for some constant ${\alpha}>-1$, for every $k$
and every ${\theta}\in\Theta$, the maximum likelihood estimator
$\hat{\theta}_{n}$ satisfies $\hat{\theta}_{n}\rightarrow\theta_{0}$ almost
surely under $\theta_{0}$.

###  4.1 Identifiability From Score Equation

For computational ease the maximum likelihood estimator may be characterized
as a solution to the likelihood equations $\dot{\iota}_{n}(f_{\theta})=0$, for
$\dot{\iota}_{n}$ given in (2.4). The asymptotic version of this function is

$\dot{\iota}(f_{\theta})=\sum_{k=1}^{\infty}\frac{\dot{f}_{\theta}}{f_{\theta}}(k)p_{>k}^{\scriptscriptstyle(0)}-\frac{\sum_{k=1}^{\infty}p_{k}^{\scriptscriptstyle(0)}\dot{f}_{\theta}(k)}{\sum_{k=1}^{\infty}p_{k}^{\scriptscriptstyle(0)}f_{\theta}(k)}.$
(4.1)

It follows from (3.8) that the true parameter ${\theta}_{0}$ solves the
equation $\dot{\iota}(f_{\theta})=0$. The following proposition shows that the
processes $\dot{\iota}_{n}$ tend uniformly to $\dot{\iota}$. The proof is
similar to the proof of Theorem 4.1 and will be omitted.

######  Proposition 4.2.

In the model $\\{f_{\theta}:\theta\in\Theta\\}$ stated in Section 1.3 with
compact parameter space $\Theta\subset\mathbb{R}^{d}$ and PA functions
satisfying $f_{\theta}(k)\leq Ck^{\beta}$ for some constants $C$ and $\beta<1$
or $f_{\theta}(k)=k+{\alpha}$ for some constant ${\alpha}>-1$, for every $k$
and every ${\theta}\in\Theta$, assume that, for some constants $C$ and
${\gamma}$,

$\displaystyle\bigl{\|}\dot{f}_{\theta}(k)\bigr{\|}$ $\displaystyle\leq
Ck\log^{\gamma}k,$
$\displaystyle\Bigl{\|}\frac{\dot{f}_{\theta}}{f_{\theta}}(k)\Bigr{\|}$
$\displaystyle\leq C\log^{\gamma}k.$

Then
$\sup_{\theta\in\Theta}|\dot{\iota}_{n}(f_{\theta})-\dot{\iota}(f_{\theta})|\rightarrow
0$, almost surely, as $n\rightarrow\infty$.

It follows that the maximum likelihood estimator is asymptotically the unique
solution to the likelihood equations in compact subsets of the parameter space
in which ${\theta}_{0}$ is the unique zero of
${\theta}\mapsto\dot{\iota}(f_{\theta})$ (compare Theorem 5.9 in van der Vaart
(2000)). However, in general proving global uniqueness turns out to be
difficult. We present the following partial results, starting with two useful
lemmas.

######  Lemma 4.3.

Suppose $(v_{k})_{k=1}^{\infty}$ is strictly decreasing with respect to $k$.
If $(p_{k})_{k=1}^{\infty}$ and $(q_{k})_{k=1}^{\infty}$ are probability
distributions on $\mathbb{N}_{+}$ such that $p_{k}\leq q_{k}$ for $k\leq K$
and $p_{k}>q_{k}$ for $k>K$, then

$\sum_{k=1}^{\infty}p_{k}v_{k}>\sum_{k=1}^{\infty}q_{k}v_{k}.$

In case $(v_{k})_{k=1}^{\infty}$ is strictly increasing, the inequality is
true in the opposite direction.

######  Lemma 4.4.

For a probability distribution $(p_{k})_{k=1}^{\infty}$ and nonnegative
sequence $(w_{k})_{k=1}^{\infty}$ such that
$\sum_{k=1}^{\infty}p_{k}w_{k}<\infty$, define $(q_{k})_{k=1}^{\infty}$ by
$q_{k}=p_{k}w_{k}/\sum_{j}p_{j}w_{j}$. If $(w_{k})_{k=1}^{\infty}$ is strictly
increasing, then there exists a $K$ such that $p_{k}\geq q_{k}$ for $k\leq K$
and $p_{k}<q_{k}$ for $k>K$. If $(w_{k})_{k=1}^{\infty}$ is strictly
decreasing, then there exists a $K$ such that $p_{k}\leq q_{k}$ for $k\leq K$
and $p_{k}>q_{k}$ for $k>K$.

We say that changing from $\theta_{0}$ to $\theta$ induces monotonicity if
$f_{\theta}(k)/f_{\theta_{0}}(k)$ is either strictly increasing or strictly
decreasing in $k$.

######  Lemma 4.5.

If changing from ${\theta}_{0}$ to ${\theta}$ induces monotonicity on
$f_{\theta}/f_{\theta_{0}}$ for every $\theta$ in a subset
$\Theta^{\prime}\subset\Theta$, then $\dot{\iota}(f_{\theta})\neq 0$, for
every $\theta\in\Theta^{\prime}$.

For illustrative purposes, we next study in detail different variants of the
parametric form $k\mapsto(k+\alpha)^{\beta}$ in the sections 4.1.1–4.1.3.
However, we point out our results work in much more general capacity, as
studied in section 4.1.5. For instance, one could easily work out that our
results apply to the parametric specification
$k\mapsto(\log(k+\alpha))^{\beta}$. We provide a toy example in Section 4.1.4
for such a parametric form besides $k\mapsto(k+\alpha)^{\beta}$ and their
variants.

####  4.1.1 The model $f_{\alpha,\beta}(k)=(k+\alpha)^{\beta}$

In the case that $f_{\alpha,\beta}(k)=(k+\alpha)^{\beta}$, changing the
parameter $(\alpha_{0},\beta_{0})$ to another parameter $({\alpha},{\beta})$
does not always induce monotonicity on
$f_{\alpha,\beta}/f_{\alpha_{0},\beta_{0}}$. An analysis of the derivative of
the function $g(x)=(x+\alpha)^{\beta}/(x+\alpha_{0})^{\beta_{0}}$ yields that
$f_{\alpha,\beta}/f_{\alpha_{0},\beta_{0}}(k)$ is increasing in $k$ on
$\\{(\alpha,\beta):\beta-\beta_{0}+\beta\alpha_{0}-\beta_{0}\alpha\geq
0,\beta\geq\beta_{0}\\}$ and is decreasing on the set
$\\{(\alpha,\beta):\beta-\beta_{0}+\beta\alpha_{0}-\beta_{0}\alpha\geq
0,\beta\leq\beta_{0}\\}$. The preceding technique does prove that
$(\alpha_{0},\beta_{0})$ is a unique root of $\dot{\iota}(f_{\alpha,\beta})=0$
in these sets, but this does not exhaust the full parameter set.

However, the lemmas 4.4 and 4.3 may be used to prove local uniqueness. The
Hessian matrix of $\iota(f_{\alpha,\beta})$ evaluated at
$(\alpha_{0},\beta_{0})$ can be calculated as, with the shorthand notation
$f_{0}:=f_{\alpha_{0},\beta_{0}}$,

$\ddot{\iota}(f_{\alpha_{0},\beta_{0}})=-\frac{1}{a^{2}}\begin{pmatrix}ab-c^{2}&ad-
ce\\\ ad-ce\ \ &af-e^{2}\end{pmatrix},$ (4.2)

where the quantities $a,b,c,d,e,f$ are defined as

$\displaystyle a$
$\displaystyle=\sum_{k=1}^{\infty}p_{k}^{\scriptscriptstyle(0)}f_{0}(k),$
$\displaystyle b$
$\displaystyle=\sum_{k=1}^{\infty}p_{k}^{\scriptscriptstyle(0)}f_{0}(k)\frac{\beta_{0}^{2}}{(k+\alpha_{0})^{2}},$
$\displaystyle c$
$\displaystyle=\sum_{k=1}^{\infty}p_{k}^{\scriptscriptstyle(0)}f_{0}(k)\frac{\beta_{0}}{k+\alpha_{0}},$
$\displaystyle d$
$\displaystyle=\sum_{k=1}^{\infty}p_{k}^{\scriptscriptstyle(0)}f_{0}(k)\frac{\beta_{0}}{k+\alpha_{0}}\log(k+\alpha_{0}),$
$\displaystyle e$
$\displaystyle=\sum_{k=1}^{\infty}p_{k}^{\scriptscriptstyle(0)}f_{0}(k)\log(k+\alpha_{0}),\qquad$
$\displaystyle f$
$\displaystyle=\sum_{k=1}^{\infty}p_{k}^{\scriptscriptstyle(0)}f_{0}(k)\log^{2}(k+\alpha_{0}).$

It follows that $ab-c^{2}$ is strictly positive, since by the Cauchy–Schwarz
inequality, where $K$ follows the law
$(p_{k}^{\scriptscriptstyle(0)})_{k=1}^{\infty}$,

$\Bigl{(}\mathbb{E}_{p_{0}}\Bigl{[}f_{0}(K)\frac{\beta_{0}}{K+\alpha_{0}}\Bigr{]}\Bigr{)}^{2}<{\mathbb{E}_{p_{0}}[f_{0}(K)]}\,{\mathbb{E}_{p_{0}}\Bigl{[}f_{0}(K)\frac{\beta_{0}^{2}}{(K+\alpha_{0})^{2}}\Bigr{]}}.$

The same arguments work to prove that $af-e^{2}>0$ and $bf-d^{2}>0$. The
determinant of the Hessian matrix is given by

$\displaystyle|\ddot{\iota}(f_{\alpha_{0},\beta_{0}})|$
$\displaystyle=a^{2}bf+c^{2}e^{2}-abe^{2}-afc^{2}-a^{2}d^{2}-c^{2}e^{2}+2acde$
$\displaystyle=a^{2}(bf-d^{2})-ae(be-cd)-ac(cf-ed).$

This can be shown to be strictly positive by showing that both $be-cd<0$ and
$cf-ed<0$. We shall prove $be<cd$; the proof that also $cf<ed$ is similar.
Define $x_{k}=p_{k}^{\scriptscriptstyle(0)}f_{0}(k)\beta_{0}/(k+\alpha_{0})$
and $u_{k}=(k+\alpha_{0})\log(k+\alpha_{0})/\beta_{0}$,
$y_{k}=x_{k}u_{k}=p_{k}^{\scriptscriptstyle(0)}f_{0}(k)\log(k+\alpha_{0})$.
Define $p_{x}(k)=x_{k}/\sum_{j}x_{j}$ and $p_{y}(k)=y_{k}/\sum_{j}y_{j}$.
Since $u_{k}$ is strictly monotone decreasing, an application of Lemma 4.4
tells us that there exists a $K$ such that $p_{x}(k)\geq p_{y}(k)$ for $k\leq
K$ and $p_{x}(k)<p_{y}(k)$ for $k>K$. Applying Lemma 4.3 with
$w_{k}=\beta_{0}/(k+\alpha_{0})$, we see that $be<cd$.

We conclude that the Hessian matrix is negative definite, so that
$(\alpha_{0},\beta_{0})$ is a unique root of $\dot{\iota}(f_{\alpha,\beta})$
in a neighborhood around $(\alpha_{0},\beta_{0})$.

####  4.1.2 Global concavity of the model
$f_{\beta}(k)=(k+\alpha_{0})^{\beta}$ with known $\alpha_{0}$

The single-parameter model $f(k)=k^{\beta}$ is a main example of sublinear PA,
treated in (Barabási, 2016, Section 5.8). We allow in addition a nonzero
offset ${\alpha}_{0}$ and consider $f_{\beta}(k)=(k+\alpha_{0})^{\beta}$, with
$\beta$ the only parameter. The limit function (4.1) reduces to

$\dot{\iota}(f_{\beta})=\sum_{k=1}^{\infty}\log(k+\alpha_{0})p_{>k}^{\scriptscriptstyle(0)}-\frac{\sum_{k=1}^{\infty}p_{k}^{\scriptscriptstyle(0)}(k+\alpha_{0})^{\beta}\log(k+\alpha_{0})}{\sum_{k=1}^{\infty}p_{k}^{\scriptscriptstyle(0)}(k+\alpha_{0})^{\beta}}.$

The second order derivative can be computed as

$\displaystyle\ddot{\iota}(f_{\beta})=-\frac{\sum_{k=1}^{\infty}p_{k}^{\scriptscriptstyle(0)}(k+\alpha_{0})^{\beta}\log^{2}(k+\alpha_{0})}{\sum_{k=1}^{\infty}p_{k}^{\scriptscriptstyle(0)}(k+\alpha_{0})^{\beta}}+\Bigl{(}\frac{\sum_{k=1}^{\infty}p_{k}^{\scriptscriptstyle(0)}(k+\alpha_{0})^{\beta}\log(k+\alpha_{0})}{\sum_{k=1}^{\infty}p_{k}^{\scriptscriptstyle(0)}(k+\alpha_{0})^{\beta}}\Bigr{)}^{2}.$

This can be seen to be strictly negative for any $\beta\in[0,1]$, as a
consequence of the Cauchy–Schwarz inequality, with $K$ following the law
$(p_{k}^{\scriptscriptstyle(0)})_{k=1}^{\infty}$,

$\displaystyle\bigl{(}\mathbb{E}_{p_{0}}\bigl{[}(K+\alpha_{0})^{\beta}\log(K+\alpha_{0})\bigr{]}\bigr{)}^{2}<\mathbb{E}_{p_{0}}\bigl{[}(K+\alpha_{0})^{\beta}\log^{2}(K+\alpha_{0})\bigr{]}\mathbb{E}_{p_{0}}(K+\alpha_{0})^{\beta}.$

Thus, the limiting log likelihood is concave and the root of the limiting
score function is unique. Another perspective is that moving the parameter
$\beta$ away from $\beta_{0}$ induces monotonicity on
$f_{\beta}/f_{\beta_{0}}(k)=(k+\alpha_{0})^{\beta-\beta_{0}}$.

####  4.1.3 Almost-global uniqueness in case of
$f(\alpha)=(k+\alpha)^{\beta_{0}}$ with known $\beta_{0}$

For any $\alpha\neq\alpha_{0}$, the function
$k\mapsto(k+\alpha)^{\beta_{0}}/(k+\alpha_{0})^{\beta_{0}}$ is monotone
increasing (when $\alpha<\alpha_{0}$) or decreasing (when
$\alpha>\alpha_{0}$). Applying Lemma 4.5, we conclude that the root is unique
in every bounded domain.

####  4.1.4 Global concavity of the model $f(\beta)=\log^{\beta}(k)$

For further illustrative purposes, we study the toy case of the PA function
being $k\mapsto\log^{\beta}(k)$, where $\beta$ is the only parameter. We
easily calculate the score function and its derivative as follows

$\displaystyle\dot{\iota}(f_{\beta})$
$\displaystyle=\sum_{k=1}^{\infty}\log\log(k){p^{\scriptscriptstyle(0)}_{>k}}-\frac{\sum_{k=1}^{\infty}p_{k}^{\scriptscriptstyle(0)}\log^{\beta}(k)\log\log(k)}{\sum_{k=1}^{\infty}p_{k}^{\scriptscriptstyle(0)}\log^{\beta}(k)},$
$\displaystyle\ddot{\iota}(f_{\beta})$
$\displaystyle=-\frac{\sum_{k=1}^{\infty}p_{k}^{\scriptscriptstyle(0)}\log^{\beta}(k)(\log\log(k))^{2}}{\sum_{k=1}^{\infty}p_{k}^{\scriptscriptstyle(0)}\log^{\beta}(k)}+\biggl{(}\frac{\sum_{k=1}^{\infty}p_{k}^{\scriptscriptstyle(0)}\log^{\beta}(k)\log\log(k)}{\sum_{k=1}^{\infty}p_{k}^{\scriptscriptstyle(0)}\log^{\beta}(k)}\biggr{)}^{2}.$

Similar to the analysis in Section 4.1.2, $\ddot{\iota}(f_{\beta})$ here is
strictly negative for any $\beta>0$ by the Cauchy–Schwarz inequality

$\displaystyle\bigl{(}\mathbb{E}_{p_{0}}\bigl{[}\log^{\beta}(K)\log\log(K)\bigr{]}\bigr{)}^{2}<\mathbb{E}_{p_{0}}\bigl{[}\log^{\beta}(K)(\log\log(K))^{2}\bigr{]}\mathbb{E}_{p_{0}}[\log^{\beta}(K)],$

where $K\sim(p_{k}^{\scriptscriptstyle(0)})_{k=1}^{\infty}$ is an auxiliary
random variable. As such, the limiting score function is monotone decreasing
with respect to the parameter $\beta$, and has a unique zero at
$\beta=\beta_{0}$. Furthermore, the limiting log-likelihood is concave with a
unique maximizer.

####  4.1.5 The case of general $f_{\theta}$

In practice, it could happen that the score function (2.4) has multiple roots,
particularly if its limit (4.1) has multiple roots. In such cases, we may
employ the empirical estimators in Gao et al. (2017) to identify the correct
one. These are defined as follows (see their Equation (2)):

$\hat{r}_{k}(n)=\frac{N_{>k}(n)}{N_{k}(n)}.$

In Gao et al. (2017) these estimators are shown to converge to
$f_{\theta_{0}}(k)/(\sum_{j=1}^{\infty}f_{\theta_{0}}(j)p^{\scriptscriptstyle(0)}_{j})$,
almost surely. This suggests an empirical estimator for the parameter
${\theta}\in\Theta\subset\mathbb{R}^{d}$ as the solution of the system of
equations:

$\frac{f_{\theta}(k)}{f_{\theta}(1)}=\frac{\hat{r}_{k}(n)}{\hat{r}_{1}(n)},\qquad
k=2,3,\ldots,d+1.$ (4.3)

Two computational strategies suggest themselves. If finding the set of
solutions to the likelihood equations $\dot{\iota}_{n}(f_{\theta})=0$ is
easier than solving (4.3), then we may select from this set the solution that
minimizes

$\sum_{k=2}^{d+1}\left|\frac{f_{\theta}(k)}{f_{\theta}(1)}-\frac{\hat{r}_{k}(n)}{\hat{r}_{1}(n)}\right|.$

On the other hand, if (4.3) is easier to solve than (2.4), then we may find
the solution of likelihood equations in a neighborhood of the solution of
(4.3), possibly by an iterative scheme such as Newton’s algorithm.

In both cases the resulting estimator will be consistent.

######  Theorem 4.6.

Under the conditions of Theorem 4.1, the solution of the likelihood equation
resulting from either of the two indicated procedures is almost surely
consistent for ${\theta}_{0}$.

##  5 Asymptotic Normality

We prove asymptotic normality of consistent solutions to the likelihood
equations, as follows. Recall that for $\theta\in\Theta\subset\mathbb{R}^{d}$
and every $k\in\mathbb{N}_{+}$, $\dot{f}_{\theta}(k)\in\mathbb{R}^{d}$ and
$\ddot{f}_{\theta}(k)\in\mathbb{R}^{d\times d}$ are the gradient and the
Hessian matrix of $f_{\theta}(k)$ with respect to $\theta$, respectively.

######  Theorem 5.1.

In the model $\\{f_{\theta}:\theta\in\Theta\\}$ stated in Section 1.3 with
compact parameter space $\Theta\subset\mathbb{R}^{d}$ and PA functions
satisfying $f_{\theta}(k)\leq Ck^{\beta}$ for some constants $C$ and $\beta<1$
or $f_{\theta}(k)=k+{\alpha}$ for some constant ${\alpha}>-1$, for every $k$,
and every ${\theta}\in\Theta$, assume in addition that the PA functions
satisfy, for some constants $\gamma>0$ and $C>0$,

$\displaystyle\bigl{\|}\dot{f}_{\theta}(k)\bigr{\|}+\bigl{\|}\ddot{f}_{\theta}(k)\bigr{\|}$
$\displaystyle\leq Ck\log^{\gamma}k,$ (5.1)
$\displaystyle\Bigl{\|}\frac{\dot{f}_{\theta}}{f_{\theta}}(k)\Bigr{\|}+\Bigl{\|}\frac{\ddot{f}_{\theta}}{f_{\theta}}(k)\Bigr{\|}$
$\displaystyle\leq C\log^{\gamma}k.$ (5.2)

Then as $n\rightarrow\infty$, for any consistent sequence of solutions
$\hat{\theta}_{n}$ to the likelihood equations
$\dot{\iota}_{n}(f_{\theta})=0$, the sequence
$\sqrt{n}(\hat{\theta}_{n}-\theta_{0})$ converges in distribution to the
$N(0,V_{0}^{-1})$ distribution, for $V_{0}$ the ($d\times d$) matrix given by

$V_{0}=\sum_{k=1}^{\infty}\frac{\dot{f}_{\theta_{0}}\dot{f}_{\theta_{0}}^{T}}{f_{\theta_{0}}^{2}}(k){p^{\scriptscriptstyle(0)}_{>k}}-\Bigl{(}\sum_{k=1}^{\infty}\frac{\dot{f}_{\theta_{0}}}{f_{\theta_{0}}}(k){p^{\scriptscriptstyle(0)}_{>k}}\Bigr{)}\Bigl{(}\sum_{k=1}^{\infty}\frac{\dot{f}_{\theta_{0}}^{T}}{f_{\theta_{0}}}(k){p^{\scriptscriptstyle(0)}_{>k}}\Bigr{)}.$
(5.3)

In particular this is true for the maximum likelihood estimator if
${\theta}_{0}$ is identifiable and interior to $\Theta$.

######  Corollary 5.2.

For the parametric family $f_{\alpha,\beta}(k)=(k+\alpha)^{\beta}$ with true
parameter $(\alpha_{0},\beta_{0})$ in the interior of the parameter set
$\Theta_{\varepsilon}=[-1+\varepsilon,\varepsilon^{-1}]\times[0,1]$ for some
small $\varepsilon>0$, the maximum likelihood estimator
$(\hat{\alpha}_{n},\hat{\beta}_{n})$ satisfies, as $n\rightarrow\infty$,

$\sqrt{n}\left(\begin{bmatrix}\hat{\alpha}_{n}\\\
\hat{\beta}_{n}\end{bmatrix}-\begin{bmatrix}\alpha_{0}\\\
\beta_{0}\end{bmatrix}\right)\rightsquigarrow
N(0,V_{\alpha_{0},\beta_{0}}^{-1}),$

where $V_{\alpha_{0},\beta_{0}}$ is defined as the negative of
$\ddot{\iota}(f_{\alpha_{0},\beta_{0}})$ in (4.2).

To test whether the PA function is affine we may compare an estimator
$\hat{\beta}_{n}$ for $\beta$ in the model considered in Corollary 5.2 to the
value ${\beta}_{0}=1$ in the affine case. The Wald-type test statistic admits
the form

$T_{n}:=\frac{{n}^{1/2}(\hat{\beta}_{n}-1)}{(V_{\hat{\alpha},1}^{-1})_{2,2}^{1/2}},$

where $V_{\hat{\alpha},1}$ is obtained by plugging in $(\hat{\alpha}_{n},1)$
in the definition of $V_{\alpha,\beta}$ in (4.2) and
$(V_{\hat{\alpha},1}^{-1})_{2,2}$ is the $(2,2)$-element of
$V_{\hat{\alpha},1}^{-1}$.

######  Theorem 5.3.

For the parametric family $f_{\alpha,\beta}(k)=(k+\alpha)^{\beta}$ with true
parameter $(\alpha_{0},1)$ in the parameter set
$\Theta_{\varepsilon}=(-1+\varepsilon,\varepsilon^{-1})\times[0,1]$ for some
small $\varepsilon>0$, the sequence $T_{n}$ tends in distribution to
$Z{1}_{Z\leq 0}$, for $Z$ a standard normal variable.

##  6 A Remedy to the History Problem

A practical problem with the maximum likelihood estimator is that the log
likelihood function (2.2) and its derivative (2.4) depend on the history of
the network evolution. In many real-world applications observing the entire
history is impossible or too costly, and only the final snapshot at time $n$
is available. For instance, when building a social network model, we may
observe the final network, but recovering how it exactly evolved into its
current shape is difficult—we would need to check with everyone when s/he
became friends with everyone else and establish a strict time order.

This problem is solvable. The log likelihood (2.2) or its derivative (2.4)
consist of two terms, and the history problem only arises in the second term,
which results from the norming constant to the likelihood. The first terms in
(2.2) or (2.4) depend on the network only through $P_{>k}(n)$ and hence are
available from the final snapshot. The second terms are the Césaro averages of
$\log S_{f_{\theta}}(t-1)$ and its derivative
$S_{\dot{f}_{\theta}}(t-1)/S_{f_{\theta}}(t-1)$, respectively. Because
$S_{h}(t)/t$ tends to the limit $\sum_{k=1}^{\infty}h(k)p_{k}$ almost surely,
as $t\rightarrow\infty$, these Césaro averages are asymptotically actually
very close to $\log S_{f_{\theta}}(n)/n$ and
$S_{\dot{f}_{\theta}}(n)/S_{f_{\theta}}(n)$, which do depend only on the
snapshot of the network at time $n$. The remedy is to replace (2.2) or (2.4)
by

$\displaystyle\tilde{\iota}_{n}(f_{\theta})$
$\displaystyle=\sum_{k=1}^{\infty}\log f_{\theta}(k)P_{>k}(n)-\log
S_{f_{\theta}}(n),$ (6.1) $\displaystyle\dot{\tilde{\iota}}_{n}(f_{\theta})$
$\displaystyle=\sum_{k=1}^{\infty}\frac{\dot{f}_{\theta}}{f_{\theta}}(k)P_{>k}(n)-\frac{S_{\dot{f}_{\theta}}(n)}{S_{f_{\theta}}(n)}.$
(6.2)

Define a pseudo maximum likelihood estimator $\tilde{\theta}_{n}$ as the
maximizer of the first function or a zero of the second. Inspection of the
proof of Theorem 4.1 readily shows that this pseudo maximum likelihood
estimator is consistent under the same conditions as the maximum likelihood
estimator.

######  Theorem 6.1.

Under the conditions of Theorem 4.1, the pseudo maximum likelihood estimator
$\tilde{\theta}_{n}$ is consistent, i.e., as $n\rightarrow\infty$,
$\tilde{\theta}_{n}\rightarrow\theta_{0}$ almost surely, under $\theta_{0}$.

We have no proof of the asymptotic normality of the pseudo maximum likelihood
estimator in the same generality as for the maximum likelihood estimator, but
we note the following general theorem, and its corollary.

Let
$F_{\theta}(k)=\bigl{(}\sum_{j=1}^{k-1}\dot{f}_{\theta}/f_{\theta}(j),\dot{f}_{\theta}(k),f_{\theta}(k)\bigr{)}^{T}$.

######  Theorem 6.2.

Assume the conditions of Theorem 4.1, and in addition assume that the sequence
of random vectors
$\sum_{k=1}^{\infty}F_{{\theta}_{0}}(k)\sqrt{n}(P_{k}(n)-p_{k}^{\scriptscriptstyle(0)})$
is asymptotically normal with mean zero and covariance $W_{0}$ and that
${\theta}_{0}$ is interior to $\Theta$. Then the sequence of pseudo maximum
likelihood estimators $\tilde{\theta}_{n}$ satisfies
$\sqrt{n}(\tilde{\theta}_{n}-{\theta}_{0})\rightsquigarrow
N(0,V_{0}^{-1}VV_{0}^{-1})$, under $\theta_{0}$ as $n\rightarrow\infty$, where
$V$ is given in the proof below and $V_{0}$ is given in (5.3).

It is plausible that the sequence of empirical degrees
$\sqrt{n}(P_{k}(n)-p_{k}^{\scriptscriptstyle(0)})$ is asymptotically normally
distributed in some generality, but this has been established only for the
affine PA function (see Móri (2002) and Resnick and Samorodnitsky (2016), or
Proposition B.2). The preceding theorem requires that certain linear
combinations of the variables over $k$ are asymptotically normal, where the
coefficients typically tend to infinity with $k$ (e.g. at the order $k\log
k$). Although this convergence requires additional bounds for large $k$, the
condition of the theorem seems plausible in general, for sublinear PA models.

In the appendix we verify the condition for the interesting case of PA
functions that are eventually constant, for which the limiting degree
distribution follows a power law with exponential cut-off (Rudas, Tóth and
Valkó (2007)). This leads to the following corollary.

######  Corollary 6.3.

In the model with PA functions $f_{\theta}$ satisfying
$f_{\theta}(k)=f_{\theta}(k\wedge K)$, for some given $K\in\mathbb{N}_{+}$ for
which the eigenvalue condition $\mathop{\rm
Re}\nolimits\lambda_{2}(A_{K})<\lambda_{1}(A_{K})/2$ (as in the appendix)
holds at $\theta_{0}$, the pseudo maximum likelihood estimator
$\tilde{\theta}_{n}$ satisfies
$\sqrt{n}(\tilde{\theta}_{n}-{\theta}_{0})\rightsquigarrow
N(0,V_{0}^{-1}VV_{0}^{-1})$, provided that it is consistent at ${\theta}_{0}$.

The limiting covariance in the preceding theorem and corollary may be
complicated. To get around this, we propose the following bootstrap procedure.
Suppose that we observe the final snapshot of the network $G_{n}$ with $n$
nodes.

  1. 1.

Obtain the pseudo maximum likelihood estimator $\tilde{\theta}_{n}$ based on
$G_{n}$.

  2. 2.

Given $\tilde{\theta}_{n}$, simulate PA networks $(G_{m}^{(i)})_{i=1}^{s}$
with PA function $f_{\tilde{\theta}_{n}}$, each with $m$ nodes.

  3. 3.

Obtain the pseudo maximum likelihood estimator $\tilde{\theta}_{m}^{(i)}$
based on $G_{m}^{(i)}$, for $i\in[s]$.

  4. 4.

Approximate the limit variance of $\sqrt{n}(\tilde{\theta}_{n}-\theta_{0})$ by

$\tilde{\Sigma}_{m,s}(f_{\tilde{\theta}_{n}}):={m}\biggl{\\{}\frac{1}{s}\sum_{i=1}^{s}\tilde{\theta}_{m}^{(i)}\bigl{(}\tilde{\theta}_{m}^{(i)}\bigr{)}^{T}-\Bigl{(}\frac{1}{s}\sum_{i=1}^{s}\tilde{\theta}_{m}^{(i)}\Bigr{)}\Bigl{(}\frac{1}{s}\sum_{i=1}^{s}\tilde{\theta}_{m}^{(i)}\Bigr{)}^{T}\biggr{\\}}.$
(6.3)

This procedure will be consistent under a mild continuity condition on the
model ${\theta}\mapsto f_{\theta}$.

Section 7 of Gao and van der Vaart (2017) on the affine PA model addressed a
different history problem—the history of the initial degrees. In the present
paper we study the case of fixed initial degree $1$, but need to know how the
incoming nodes connect, a problem that did not arise in the affine model. The
empirical estimator considered in Gao et al. (2017) is curiously free of the
history problem.

##  7 Connecting the empirical estimator and the pseudo maximum likelihood
estimator

Write $f(k)$ as $\theta_{k}$ and suppose that we are interested in estimating
the _infinite-dimensional_ vector $(\theta_{k})_{k=1}^{\infty}$. The pseudo
log-likelihood function (6.1) then takes the form

$\tilde{\iota}_{n}(\theta)=\sum_{k=1}^{\infty}(\log\theta_{k})P_{>k}(n)-\log
S_{\theta}(n),$

where $S_{\theta}(t)=\sum_{k=1}^{\infty}\theta_{k}N_{k}(t)$ is the total
preference. Taking the derivative with respect to $\theta_{k}$, we obtain

$\frac{\partial\tilde{\iota}_{n}(\theta)}{\partial\theta_{k}}=\frac{P_{>k}(n)}{\theta_{k}}-\frac{N_{k}(n)}{\sum_{j=1}^{\infty}\theta_{j}N_{j}(n)}.$

Setting this equation to zero for every $k$, we obtain the system of equations
(for simplicity assume that $N_{k}(n)>0$ for every $k$)

$\frac{\theta_{k}}{\sum_{j=1}^{\infty}\theta_{j}P_{j}(n)}=\frac{N_{>k}(n)}{N_{k}(n)},\qquad
k\in\mathbb{N}_{+}.$

The right side is the aforementioned empirical estimator $\hat{r}_{k}(n)$,
defined in Gao et al. (2017). The PA function $f$ or parameter ${\theta}_{k}$
is identifiable up to a scale factor only. We conclude that the empirical
estimator is the pseudo maximum likelihood estimator if we do not impose any
parametric assumption and wish to estimate $\theta_{k}$ individually for any
$k$.

##  8 Numerical Illustrations

In this section we numerically study the performance of the maximum likelihood
estimator, the pseudo maximum likelihood estimator, the empirical estimator,
the Wald test and the bootstrap estimator. We simulated data using the
following three examples of a PA function, each of the type
$f(k)=(k+{\alpha})^{\beta}$:

$\displaystyle f^{\scriptscriptstyle(2)}(k)$ $\displaystyle=k^{2/3},$
$\displaystyle{\alpha}=0,$ $\displaystyle\ {\beta}=2/3,$ $\displaystyle
f^{\scriptscriptstyle(4)}(k)$ $\displaystyle=(k+4)^{4/5},$
$\displaystyle{\alpha}=4,$ $\displaystyle\ {\beta}=4/5,$ $\displaystyle
f^{\scriptscriptstyle(5)}(k)$ $\displaystyle=k+2,$ $\displaystyle{\alpha}=2,$
$\displaystyle\ \beta=1.$

In every setting we conducted $N=1000$ repetitions of the experiment in which
we simulated a PA tree of $n$ nodes, for varying $n$, and computed the
estimators and/or test. For every estimator $\hat{\theta}_{i}$ of
$\theta=({\alpha},{\beta})$ we computed the sample mean difference
$(1/N)\sum_{i=1}^{N}(\hat{\theta}_{i}-\theta_{0})$ and the rescaled sample
covariance matrix
$(n/N)\sum_{i=1}^{N}(\hat{\theta}_{i}-\theta_{0})(\hat{\theta}_{i}-\theta_{0})^{T}$.
According to our theory, for large $n$ the first should be close to zero, and
the second should be close to the deterministic matrix $V_{0}^{-1}$, for
$V_{0}$ given in (5.3), which depends on the PA function.

###  8.1 MLE

The limiting covariance matrix $V_{0}^{-1}$ of
$\sqrt{n}(\hat{\theta}_{n}-\theta_{0})$ under the PA function
$f^{\scriptscriptstyle(2)}$ is computed to be

$\begin{pmatrix}169.30\ &47.56\\\ 47.56&14.94\end{pmatrix}.$

Table 1 gives the results of our simulation experiments, where we simulated PA
trees of three different sizes: $n=10^{4}$, $n=10^{5}$ or $n=10^{6}$ nodes. As
expected the sample mean difference decreases with $n$, while the sample
covariance swings around the expected limit.

# of nodes | Sample Mean Difference | Sample Covariance
---|---|---
$10^{4}$ |  | 7.29e-03 | 3.24e-05
---|---
| 173.19 | 48.14
---|---
48.14 | 15.08
$10^{5}$ |  | 5.71e-04 | -7.03e-05
---|---
| 167.89 | 47.23
---|---
47.23 | 15.03
$10^{6}$ |  | 1.71e-04 | -2.59e-05
---|---
| 163.02 | 46.60
---|---
46.60 | 14.97
Table 1: Sample mean difference and covariance of the maximum likelihood
estimator of the parameters ${\alpha}$ and ${\beta}$ of the PA function
$f(k)=(k+{\alpha})^{\beta}$ for trees of three sizes of node sets $n$
generated according to PA function $f^{\scriptscriptstyle(2)}$.

We conducted the same experiment with PA functions $f^{\scriptscriptstyle(4)}$
and $f^{\scriptscriptstyle(5)}$, but for trees of a single size of $n=10^{6}$
nodes. The results presented in Table 2 again confirm the theory.

PA Function | Sample Mean Difference | Sample Covariance | Limit Variance
---|---|---|---
$f^{\scriptscriptstyle(4)}$ |  | -2.43e-02 | -2.72e-03
---|---
| 42764.75 | 4743.46
---|---
4743.46 | 540.71
| 42429.33 | 4716.76
---|---
4716.76 | 539.75
$f^{\scriptscriptstyle(5)}$ |  | -5.35e-04 | -1.01e-04
---|---
| 1817.94 | 325.66
---|---
325.66 | 63.07
| 1762.05 | 316.58
---|---
316.58 | 61.64
Table 2: Sample mean difference and covariance of the maximum likelihood
estimator of the parameters ${\alpha}$ and ${\beta}$ of the PA function
$f(k)=(k+{\alpha})^{\beta}$ for trees of $n=10^{6}$ nodes generated according
to the PA functions $f^{\scriptscriptstyle(4)}$ and
$f^{\scriptscriptstyle(5)}$. The last column gives the covariance matrix
$V_{0}^{-1}$.

###  8.2 Comparing the estimators

We compared the performance of the empirical estimator (4.3), the pseudo
maximum likelihood estimator and the maximum likelihood estimator on samples
of trees of $n=10^{6}$ nodes generated using the PA function
$f^{\scriptscriptstyle(2)}$. Table 3 gives numerical results, while Figure 1
shows QQ-plots of the estimators. The first line of Table 3 repeats the
relevant (third) line of Table 1.

The comparison is particularly interesting as the theoretical variances of the
empirical estimator and the pseudo maximum likelihood estimator are unknown.
In our experiment the sample covariance of the pseudo maximum likelihood
estimator is only twice bigger than that of the maximum likelihood estimator.
In contrast, the sample covariance of the empirical estimator is larger than
that of the other two estimators by an order of magnitude. We conclude that it
helps to use a parametric model, if this can be correctly specified.

As is evident from Figure 1, the maximum likelihood estimator and pseudo
maximum likelihood estimator are asymptotically normal. The same seems to be
true for the empirical estimator, with the largest visible possible deviation
in the left tail of its distribution.

Estimator | Sample Mean Difference | Sample Covariance
---|---|---
MLE |  | 1.71e-04 | -2.59e-05
---|---
| 163.02 | 46.60
---|---
46.60 | 14.97
EE |  | 1.02e-03 | 2.58e-04
---|---
| 6840.42 | 2965.18
---|---
2965.18 | 1297.15
PMLE |  | -1.05e-03 | -3.79e-04
---|---
| 297.11 | 85.40
---|---
85.40 | 26.20
Table 3: Comparison of the maximum likelihood estimator, the empirical
estimator and the pseudo maximum likelihood estimator for estimating the
parameters ${\alpha}$ and ${\beta}$ of the PA function
$f(k)=(k+{\alpha})^{\beta}$ based on trees with $n=10^{6}$ nodes generated
according to PA function $f^{\scriptscriptstyle(2)}$.

(a) Estimators of $\alpha$

(b) Estimators of $\beta$

Figure 1: QQ-plots of 1000 realizations of the maximum likelihood estimator,
the empirical estimator and the pseudo maximum likelihood estimator of the
parameters ${\alpha}$ and ${\beta}$ of the PA function
$f(k)=(k+{\alpha})^{\beta}$ based on trees with $n=10^{6}$ nodes generated
according to PA function $f^{\scriptscriptstyle(2)}$.

###  8.3 The Wald-type test for affinity

We applied the Wald-type test for affinity given in Corollary 5.3 on PA trees
generated using $f^{\scriptscriptstyle(4)}$ and $f^{\scriptscriptstyle(5)}$.
The nominal size of the tests was set at 0.05. In both cases we registered the
proportion of repetitions in which the null hypothesis was rejected. As
$f^{\scriptscriptstyle(5)}$ is affine, rejection constitutes a type-I error in
this case, while for $f^{\scriptscriptstyle(4)}$ the proportion of rejections
is a measure of the test’s power. The results are summarized in Table 4. The
expected proportion of type-I errors for $f^{\scriptscriptstyle(5)}$ was 5.2%,
close to the nominal value. The power of the test at
$f^{\scriptscriptstyle(4)}$ was outright 1, resulting from the fact that the
number of nodes was large and $f^{\scriptscriptstyle(4)}$ is far from affine.

PA function | Type-I error | Power
---|---|---
$f^{\scriptscriptstyle(4)}$ | - | 1
$f^{\scriptscriptstyle(5)}$ | .052 | -
Table 4: Proportion of rejections of the Wald-type affinity test with nominal
size 0.05 based on PA trees with $n=10^{6}$ nodes for two PA functions.

###  8.4 The pseudo maximum likelihood estimator with bootstrapped variance

We conducted two experiments to illuminate the bootstrap procedure for the
pseudo maximum likelihood estimator in Section 6, applied to the model
$f(k)=(k+{\alpha})^{\beta}$. In both cases the data was generated under the PA
function $f^{\scriptscriptstyle(2)}$, corresponding to ${\alpha}=0$ and
${\beta}=2/3$. Furthermore, the bootstrap sample size in (6.3) was set to
$m=10^{5}$ and the number of bootstrap replicates to $s=10^{3}$.

In the first experiment we performed Wald type tests for the hypotheses
$H_{0}:{\alpha}=0$ and $H_{0}:{\beta}=2/3$ using the relevant coordinate of
the pseudo maximum likelihood estimator
$(\tilde{\alpha}_{n},\tilde{\beta}_{n})$ with its variance estimated by the
bootstrap estimator (6.3). In the case of the hypothesis $H_{0}:{\alpha}=0$
this entails the test statistic

$T_{n,m,s}:=\frac{n^{1/2}\,\tilde{\alpha}_{n}}{\bigl{(}\tilde{\Sigma}_{m,s}(f_{\alpha_{0},\tilde{\beta}_{n}})_{1,1}\bigr{)}^{1/2}}.$

The test for $H_{0}:{\beta}=2/3$ is similarly based on a standardized version
of $\tilde{\beta}_{n}$. We rejected the null hypothesis when
$|T_{n,m,s}|>z_{0.025}$, corresponding to the working hypothesis that
$T_{n,m,s}\sim N(0,1)$ under the null hypothesis and nominal size $0.05$.
Table 5 gives the proportions of rejections in $N=1000$ repetitions of the
experiment.

$H_{0}$ | Type-I error
---|---
$\alpha=0$ | 0.047
$\beta=2/3$ | 0.063
Table 5: Proportions of rejections of the size-0.05 Wald-type tests with
bootstrapped variance of the null hypotheses $H_{0}:{\alpha}=0$ and
$H_{0}:{\beta}=2/3$ on the PA function $f(k)=(k+{\alpha})^{\beta}$ for trees
with $n=10^{6}$ nodes generated according to PA function
$f^{\scriptscriptstyle(2)}$. The bootstrap sample size was set to $m=10^{5}$
and the number of bootstrap replicates to $s=10^{3}$.

In the second experiment we simulated $N=1000$ replicates of the normalized
and projected pseudo maximum likelihood estimator

$d_{n}:=r_{n}^{T}\sqrt{n}\bigl{(}\tilde{\Sigma}_{m,s}(f_{\tilde{\alpha}_{n},\tilde{\beta}_{n}})\bigr{)}^{-1/2}\binom{\tilde{\alpha}_{n}}{\tilde{\beta}_{n}},$

where $r_{n}$ are i.i.d. random vectors from the unit circle. Figure 2 shows a
QQ-plot of these 1000 values against the standard normal distribution.

Figure 2: QQ-plot of 1000 realizations of the randomly projected pseudo
maximum likelihood estimators in the model $f(k)=(k+{\alpha})^{\beta}$
normalized by their bootstrapped variances based on trees with $n=10^{6}$
nodes generated according to $f^{\scriptscriptstyle(2)}$. The bootstrap sample
size was set to $m=10^{5}$ and the number of bootstrap replicates to
$s=10^{3}$.

Both experiments support the practical use of asymptotic normal limit theory
for the pseudo maximum likelihood estimator and validate the bootstrap
procedure for estimating the asymptotic variance.

##  9 Proofs of main results

###### Proof of Proposition 3.1.

We apply the general results due to Jagers (1975) and Nerman (1981), as
summarized in the appendix, Section A.

The events of the pure birth process (3.1) can be represented as
$T_{1}<T_{1}+T_{2}<T_{1}+T_{2}+T_{3}<\cdots$, for independent exponential
random variables $(T_{k})_{k=1}^{\infty}$ with rates $(f(k))_{k=1}^{\infty}$.
The total number of births $\xi(t)=\int{1}_{(0,t]}(u)\,\xi(du)$ at time $t$ is
equal to $\sum_{l=1}^{\infty}{1}_{(0,t]}(T_{1}+\cdots+T_{l})$, which tends to
infinity almost surely as $t\rightarrow\infty$ by the assumption that the PA
function is bounded from below by $f(1)>0$. The birth times are clearly not
restricted to any lattice, and the functions
$t\mapsto\mathbb{E}[\varphi_{i}(t)]$ are continuous almost everywhere in view
of their monotonicity. Thus, it suffices to show that the Malthusian parameter
exists and to verify conditions (A.9) and (A.10).

Since $\int\mathrm{e}^{-\lambda
u}\,\xi(du)=\sum_{l=1}^{\infty}\mathrm{e}^{-\lambda(T_{1}+\cdots+T_{l})}$,

$\displaystyle\mathbb{E}\Bigl{[}\int_{0}^{\infty}\mathrm{e}^{-\lambda
u}\,\xi(du)\Bigr{]}$
$\displaystyle=\mathbb{E}\Bigl{[}\sum_{l=1}^{\infty}\mathrm{e}^{-\lambda(T_{1}+\cdots+T_{l})}\Bigr{]}=\sum_{l=1}^{\infty}\prod_{i=1}^{l}\frac{f(i)}{\lambda+f(i)},$
$\displaystyle\mathbb{E}\Bigl{[}\int_{0}^{\infty}\mathrm{e}^{-\lambda
u}\xi(du)\Bigr{]}^{2}$
$\displaystyle=\mathbb{E}\Bigl{[}\sum_{l=1}^{\infty}\mathrm{e}^{-\lambda(T_{1}+\dots+T_{l})}\Bigr{]}^{2}=\sum_{k=1}^{\infty}\sum_{l=1}^{\infty}\mathbb{E}\bigl{[}\mathrm{e}^{-\lambda(T_{1}+\cdots+T_{k})-\lambda(T_{1}+\cdots+T_{l})}\bigr{]}$
$\displaystyle=\rho_{f}(2\lambda)+2\sum_{k=1}^{\infty}\sum_{l=k+1}^{\infty}\prod_{i=1}^{k}\frac{f(i)}{2\lambda+f(i)}\prod_{i=k+1}^{l}\frac{f(i)}{\lambda+f(i)}.$

The left side of the first formula is $\int_{0}^{\infty}\mathrm{e}^{-\lambda
u}\,\mu(du)=\rho_{f}(\lambda)$, and hence this formula verifies equation
(3.2).

First consider the strictly sublinear case, where $f(k)\leq Ck^{\beta}$.
Because both expressions in the display are monotone in $f$, we obtain upper
bounds on these expressions by evaluating their right sides for
$f(k)=Ck^{\beta}$. Because $\log(1+x)\geq x/2$ for $x\in[0,1]$, we have
$f(i)/(\lambda+f(i))\leq\exp(-\lambda/(2f(i)))$, for $i\geq
I_{\lambda}:=(\lambda/C)^{1/\beta}$. Since the quotients and their products
are also bounded by 1,

$\sum_{l=1}^{\infty}\prod_{i=1}^{l}\frac{f(i)}{\lambda+f(i)}\leq
I_{\lambda}+\sum_{l>I_{\lambda}}e^{-\sum_{I_{\lambda}<i\leq
l}\lambda/(2Ci^{\beta})}.$

Because $\beta<1$, the sum in the exponent is of the order $l^{1-\beta}$ as
$l\rightarrow\infty$, and hence the series is finite for every $\lambda>0$.
The function $\rho_{f}$ is monotone decreasing and tends to zero as
$\lambda\rightarrow\infty$, by the dominated convergence theorem, and to
infinity as $\lambda\downarrow 0$, by the monotone convergence theorem. It
follows that the Malthusian parameter ${\lambda^{*}}$ exists and is contained
in $(0,\infty)$. We also conclude that (A.9) is satisfied, for any
$\lambda\in(0,{\lambda^{*}})$.

The second moment of $\int_{0}^{\infty}\mathrm{e}^{-\lambda u}\,\xi(du)$ is
also finite for every $\lambda>0$. To see this, it suffices to bound the
double sum on the right. By the same arguments, for $k>I_{\lambda}$,

$\displaystyle\sum_{l=k+1}^{\infty}\prod_{i=k+1}^{l}\frac{f(i)}{\lambda+f(i)}$
$\displaystyle\leq\sum_{l=k+1}^{\infty}\mathrm{e}^{-\sum_{i=k+1}^{l}\lambda/(2Ci^{\beta})}\leq\sum_{l=k+1}^{\infty}\mathrm{e}^{-(\lambda/2C(1-\beta))[l^{1-\beta}-(k+1)^{1-\beta}]}.$

This is bounded from above by a constant for every $k$ (in fact by a multiple
of $k^{-\beta}$). Inserting this bound in the double sum, we are left with a
single sum of the same type as before, except that $\lambda$ is replaced by
$2\lambda$. We conclude that the second moment of
$\int_{0}^{\infty}\mathrm{e}^{-\lambda u}\xi(du)$ is finite for every
$\lambda>0$ as well.

We have $\mathrm{e}^{-\lambda t}\xi(t)=\int_{0}^{t}\mathrm{e}^{-\lambda
t}\xi(du)\leq\int_{0}^{t}\mathrm{e}^{-\lambda u}\,\xi(du)$, for every $t$ and
$\lambda>0$. Combined with the assumption $\varphi_{i}(t)\leq C\xi(t)^{2}$, we
see that $\mathrm{e}^{-\lambda t}\varphi_{i}(t)\leq
C\bigl{[}\int_{0}^{t}\mathrm{e}^{-\lambda u/2}\,\xi(du)\bigr{]}^{2}$. It
follows that $\mathbb{E}\sup_{t>0}\bigl{(}\mathrm{e}^{-\lambda
t}\varphi_{i}(t)\bigr{)}\leq
C\mathbb{E}\bigl{[}\int_{0}^{\infty}\mathrm{e}^{-\lambda
u/2}\,\xi(du)\bigr{]}^{2}<\infty$, for every $\lambda>0$, thus verifying
(A.10).

In the case that $f(k)=k+\alpha$, for $\alpha>-1$, the function $\rho_{f}$
takes the form $\rho_{f}(\lambda)=(1+\alpha)/(\lambda-1)$ (see Rudas, Tóth and
Valkó (2007)) and hence the Malthusian parameter is ${\lambda^{*}}=2+\alpha$
and the function $\rho_{f}$ is finite for $\lambda>1$. We conclude that (A.9)
is satisfied. We show below that the second moment of
$\int_{0}^{\infty}\mathrm{e}^{-\lambda u}\xi(du)$ is finite for every
$\lambda>1$ as well. If $\varphi_{i}(t)\leq C\xi(t)^{r}$, then
$\sup_{t>0}\mathrm{e}^{-\lambda t}\varphi_{i}(t)\leq
C\bigl{[}\int_{0}^{\infty}\mathrm{e}^{-\lambda u/r}\,\xi(du)\bigr{]}^{r}$,
which then has a finite moment for any $\lambda$ such that $\lambda/r>1$ for
some $r\leq 2$. Since $r$ can be chosen such that $r\leq 2$ and $r<2+\alpha$,
where $2+\alpha>1$, there exists ${\lambda}\in(1,{\lambda^{*}})$ that
satisfies these restrictions, and hence (A.10) follows.

For $i\rightarrow\infty$, we have
$f(i)/({\lambda}+f(i))=\exp(-\log(1+\lambda/(i+\alpha))\asymp\exp(-{\lambda}/(i+\alpha))$.
Hence, for $k\rightarrow\infty$, and ${\lambda}>1$,

$\displaystyle\sum_{l=k+1}^{\infty}\prod_{i=k+1}^{l}\frac{f(i)}{\lambda+f(i)}$
$\displaystyle\asymp\sum_{l=k+1}^{\infty}\mathrm{e}^{-\lambda\log\frac{l+\alpha}{k+\alpha}}=\sum_{l=k+1}^{\infty}\left(\frac{k+\alpha}{l+\alpha}\right)^{\lambda}\asymp
k,$ $\displaystyle\prod_{i=1}^{k}\frac{f(i)}{2\lambda+f(i)}$
$\displaystyle\asymp\mathrm{e}^{-2\lambda\log(k+{\alpha})}\asymp
k^{-2\lambda}.$

Therefore, the double sum in the second moment of
$\int_{0}^{\infty}\mathrm{e}^{-\lambda u}\xi(du)$ is bounded by a multiple of
$\sum_{k=1}^{\infty}k^{-(2\lambda-1)}$, which is finite for $\lambda>1$.

###### Proof of Corollary 3.2.

The random characteristics $\varphi_{1}(t)=h(\xi(t)+1)$ and
$\varphi_{2}(t)={1}_{\\{t\geq 0\\}}$ satisfy the conditions of Proposition
3.1, and hence $Z_{t}^{\varphi_{1}}/Z_{t}^{\varphi_{2}}$ converges almost
surely to the limit given in the proposition, as $t\rightarrow\infty$. If
$\tau_{n}$ is the first time that the PA tree $\Upsilon_{t}$ possesses $n$
nodes, then $\tau_{n}\rightarrow\infty$ almost surely, as the total number of
individuals at any given $t$ is finite, almost surely. (At any given (finite)
time, every individual has a finite number of offspring and there can be at
most finitely many generations, since a new generation can be formed not
faster than an exponential waiting time with mean $1/f(1)$.) Hence, the
sequence $Z_{\tau_{n}}^{\varphi_{1}}/Z_{\tau_{n}}^{\varphi_{2}}$ converges to
the same limit, as $n\rightarrow\infty$.

The process $Z_{t}^{\varphi_{2}}$ simply counts the number of nodes $v$ in the
tree at time $t$, and hence $Z_{\tau_{n}}^{\varphi_{2}}=n$. Since
$\xi_{v}(t)+1$ is the degree of node $v$ at time $t$, the quotient can be
rewritten as

$\frac{Z^{\varphi_{1}}_{\tau_{n}}}{Z^{\varphi_{2}}_{\tau_{n}}}=\frac{1}{n}{\sum_{v\in\Upsilon_{\tau_{n}}}h\bigl{(}\deg(v,\Upsilon_{\tau_{n}})\bigr{)}}=\frac{1}{n}{\sum_{k=1}^{\infty}h(k)N_{k}(n)}=\sum_{k=1}^{\infty}h(k)P_{k}(n).$

This is the left side of the corollary, and it remains to identify the limit
given in Proposition 3.1 as its right side. The latter is equal to

$\displaystyle\frac{\int_{0}^{\infty}\mathrm{e}^{-{\lambda^{*}}t}\mathbb{E}h(\xi(t)+1)\,dt}{\int_{0}^{\infty}\mathrm{e}^{-{\lambda^{*}}t}\mathbb{E}{1}_{\\{t\geq
0\\}}\,dt}=\frac{\sum_{k=1}^{\infty}h(k){\lambda^{*}}\int_{0}^{\infty}\mathrm{e}^{-{\lambda^{*}}t}\mathbb{P}(\xi(t)+1=k)\,dt}{{\lambda^{*}}\int_{0}^{\infty}\mathrm{e}^{-{\lambda^{*}}t}\,dt}.$

The denominator is simply $1$ and the numerator is
$\sum_{k=1}^{\infty}h(k)p_{k}$ by identifying $p_{k}$ from Equation (17) in
Gao et al. (2017).

###### Proof of Theorem 4.1.

Recall the notation $\iota_{n}(f_{\theta})$ for the scaled log likelihood in
(2.2). We show below that
$\sup_{\theta}\bigl{|}\iota_{n}(f_{\theta})-\iota(f_{\theta})-c_{n}\bigr{|}\rightarrow
0$, almost surely, for $c_{n}=n^{-1}\sum_{t=2}^{n}\log(t-1)$ and

$\iota(f_{\theta})=\sum_{k=1}^{\infty}\bigl{(}\log
f_{\theta}(k)\bigr{)}\,p_{>k}^{\scriptscriptstyle(0)}-\log\Bigl{(}\sum_{k=1}^{\infty}f_{\theta}(k)p_{k}^{\scriptscriptstyle(0)}\Bigr{)}.$

Because $\hat{\theta}_{n}$ maximizes ${\theta}\mapsto\iota_{n}(f_{\theta})$,
it next suffices to show that the limit function
${\theta}\mapsto\iota(f_{\theta})$ possesses ${\theta}_{0}$ as a well-
separated point of maximum, in the sense of van der Vaart (2000), Theorem 5.7.
Since the function ${\theta}\mapsto\iota(f_{\theta})$ is continuous on the
compact set $\Theta$, this is equivalent to showing that ${\theta}_{0}$ is a
unique point of maximum.

Consider the probability distributions
$q^{\theta}=(q_{k}^{\theta})_{k=1}^{\infty}$ on $\mathbb{N}_{+}$ defined by

$q_{k}^{\theta}=\frac{1}{c({\theta})}\frac{f_{\theta}(k)}{f_{{\theta}_{0}}(k)}p_{>k}^{\scriptscriptstyle(0)},$
(9.1)

where the norming constant satisfies

$c({\theta})=\sum_{k}\frac{f_{\theta}(k)}{f_{{\theta}_{0}}(k)}p_{>k}^{\scriptscriptstyle(0)}=\frac{\sum_{k}f_{\theta}(k)p_{k}^{\scriptscriptstyle(0)}}{\sum_{k}f_{{\theta}_{0}}(k)p_{k}^{\scriptscriptstyle(0)}},$

by (3.8). In particular $c({\theta}_{0})=1$ and hence
$q^{{\theta}_{0}}_{k}=p_{>k}^{\scriptscriptstyle(0)}$. The Kullback–Leibler
divergence of $q^{\theta}$ relative to $q^{{\theta}_{0}}$ can be seen to be
equal to $\iota(f_{{\theta}_{0}})-\iota(f_{\theta})$, and is strictly positive
unless $q^{\theta}=q^{{\theta}_{0}}$. The latter is equivalent to
$f_{\theta}\propto f_{{\theta}_{0}}$, which is excluded by the identifiability
assumption.

We finish by proving the uniform convergence, where we show that the two
components of $\iota_{n}(f_{\theta})$ converge to the two components of
$\iota(f_{\theta})$. By Fubini’s theorem the difference of the first
components satisfies

$\displaystyle\Bigl{|}\sum_{k=1}^{\infty}\log
f_{\theta}(k)(P_{>k}(n)-p_{>k}^{\scriptscriptstyle(0)})\Bigr{|}$
$\displaystyle=\Bigl{|}\sum_{j=1}^{\infty}\sum_{k=1}^{j-1}\log
f_{\theta}(k)(P_{j}(n)-p_{j}^{\scriptscriptstyle(0)})\Bigr{|}$
$\displaystyle\leq\sum_{j=1}^{\infty}\Bigl{|}\sum_{k=1}^{j-1}\log
f_{\theta}(k)\Bigr{|}\,|p_{j}^{\scriptscriptstyle(0)}-P_{j}(n)|.$

Because $\log f_{\theta}(k)\leq C\log k$ and $|u|=u+2u_{-}$, for every
$u\in\mathbb{R}$, this is bounded above by a multiple of

$\sum_{j=1}^{\infty}j\,\log
j(p_{j}^{\scriptscriptstyle(0)}-P_{j}(n))+2\sum_{j=1}^{\infty}j\log
j\,(p_{j}^{\scriptscriptstyle(0)}-P_{j}(n))_{-}.$

The first series on the left tends to zero almost surely by Corollary 3.2
applied with the function $h:j\mapsto j\log j$. The terms of the second series
are bounded from above by $j\log j\,p_{j}^{\scriptscriptstyle(0)}$ and tend to
zero almost surely for every $j$ as $n\to\infty$ by Corollary 3.2 applied with
$h:k\mapsto{1}_{\\{k\\}}$. Furthermore, the series $\sum_{j}j\log
j\,p_{j}^{\scriptscriptstyle(0)}$ converges. Hence, the second series tends to
zero almost surely by the dominated convergence theorem. The two series give
an upper bound independent of ${\theta}$ and hence the supremum over
${\theta}$ of the left side of the second last equation tends to zero almost
surely.

For the second component of $\iota_{n}(f_{\theta})$, we first note that for
some constant $C$ not depending on $\theta$,

$\displaystyle\Bigl{|}\frac{S_{f_{\theta}}(n)}{n}-\sum_{k=1}^{\infty}f_{\theta}(k)p_{k}^{\scriptscriptstyle(0)}\Bigr{|}$
$\displaystyle=\Bigl{|}\sum_{k=1}^{\infty}f_{\theta}(k)(P_{k}(n)-p_{k}^{\scriptscriptstyle(0)})\Bigr{|}\leq
C\sum_{k=1}^{\infty}k\,|p_{k}^{\scriptscriptstyle(0)}-P_{k}(n)|.$

The right side tends to zero almost surely as $n\rightarrow\infty$, as in
preceding paragraph. It follows that the supremum over ${\theta}$ of the left
side tends to zero almost surely. The limit
$\sum_{k}f_{\theta}(k)p_{k}^{\scriptscriptstyle(0)}$ is a continuous, positive
function and hence is bounded away from 0. By the continuous mapping theorem,

$\sup_{\theta}\Bigl{|}\log\frac{S_{f_{\theta}}(n)}{n}-\log\Bigl{(}\sum_{k=1}^{\infty}f_{\theta}(k)p_{k}^{\scriptscriptstyle(0)}\Bigr{)}\Bigr{|}\xrightarrow{\text{a.s.}}0.$

By Lemma 9.1 the Césaro’s sums
$n^{-1}\sum_{t=2}^{n}\log(S_{f_{\theta}}(t-1)/(t-1))=n^{-1}\sum_{t=2}^{n}\log
S_{f_{\theta}}(t-1)-c_{n}$ have the same limit.

######  Lemma 9.1 (Uniform Cesàro convergence for processes).

If $Z_{\theta}(1),Z_{\theta}(2),\ldots$ are bounded stochastic processes such
that $\sup_{\theta}|Z_{\theta}(n)-Z_{\theta}|\rightarrow 0$ almost surely, for
some process $Z_{\theta}$, then the Cesàro means
$\bar{Z}_{\theta}(n)=n^{-1}\sum_{t=1}^{n}Z_{\theta}(t)$ of $Z_{\theta}(n)$
satisfy $\sup_{\theta}|\bar{Z}_{\theta}(n)-Z_{\theta}|\rightarrow 0$, almost
surely.

###### Proof of Lemma 9.1.

For every $m$ the difference $|\bar{Z}_{\theta}(n)-Z_{\theta}|$ is bounded
from above by

$\frac{m}{n}\max_{1\leq t\leq
m}|Z_{\theta}(t)-Z_{\theta}|+\frac{1}{n}\sum_{t=m+1}^{n}|Z_{\theta}(t)-Z_{\theta}|.$

For every ${\varepsilon}>0$ there exists $m$ so that every term in the second
sum is bounded from above by ${\varepsilon}$ and hence this sum divided by $n$
is bounded by ${\varepsilon}$. The first term tends to zero as
$n\rightarrow\infty$, for every fixed $m$. This argument is true also after
taking the supremum over ${\theta}$ across.

###### Proof of Lemma 4.3.

Since $\sum_{k=1}^{\infty}p_{k}=1=\sum_{k=1}^{\infty}q_{k}$, we have
$\sum_{k=1}^{K}(q_{k}-p_{k})=\sum_{k=K+1}^{\infty}(p_{k}-q_{k})$, where the
terms of the sums are nonnegative by assumption. The strict monotonicity of
$v_{k}$ gives
$\sum_{k=1}^{K}(q_{k}-p_{k})v_{k}<\sum_{k=K}^{\infty}(p_{k}-q_{k})v_{k}$.
Rearranging the terms gives the desired result.

###### Proof of Lemma 4.4.

Since both sequences sum to $1$, it is impossible that $p_{k}>q_{k}$ for every
$k\in\mathbb{N}_{+}$. If $(w_{k})_{k=1}^{\infty}$ is strictly increasing and
$p_{j}\leq q_{j}$, then

$\frac{q_{j+1}}{p_{j+1}}=\frac{w_{j+1}}{\sum_{k=1}^{\infty}p_{k}w_{k}}>\frac{w_{j}}{\sum_{k=1}^{\infty}p_{k}w_{k}}=\frac{q_{j}}{p_{j}}\geq
1.$

By mathematical induction, $p_{k}<q_{k}$ for every $k>j$. Then $p_{k}<q_{k}$
for any $k>K$ and $K+1$ the smallest value $j$ with $p_{j}<q_{j}$ (which
cannot be $j=1$).

In the case that $(w_{k})_{k=1}^{\infty}$ is strictly decreasing, the sequence
$w^{-1}_{k}$ is strictly increasing, and we apply the preceding argument with
$p_{k}=q_{k}w_{k}^{-1}/\sum_{j}q_{j}w_{j}^{-1}$ and the roles of $p_{k}$ and
$q_{k}$ swapped.

###### Proof of Lemma 4.5.

A more illustrative view of (4.1) is as follows:

$\displaystyle\dot{\iota}(f_{\theta})$
$\displaystyle=\sum_{k=1}^{\infty}\frac{\dot{f}_{\theta}}{f_{\theta}}(k)p_{>k}^{\scriptscriptstyle(0)}-\sum_{k=1}^{\infty}\frac{p_{>k}^{0}f_{\theta}(k)/f_{\theta_{0}}(k)}{\sum_{j=1}^{\infty}{p_{>j}^{\scriptscriptstyle(0)}f_{\theta}(j)/f_{\theta_{0}}(j)}}\frac{\dot{f}_{\theta}}{f_{\theta}}(k)$
(9.2)
$\displaystyle=\sum_{k=1}^{\infty}p_{>k}^{\scriptscriptstyle(0)}\frac{\dot{f}_{\theta}}{f_{\theta}}(k)-\sum_{k=1}^{\infty}q_{k}^{\scriptscriptstyle(0,\theta)}\frac{\dot{f}_{\theta}}{f_{\theta}}(k),$

where $q^{\scriptscriptstyle(0,\theta)}_{k}\propto
p_{>k}^{\scriptscriptstyle(0)}f_{\theta}(k)/f_{\theta_{0}}(k)$ is the
probability distribution generated by reweighting
$(p_{>k}^{\scriptscriptstyle(0)})_{k=1}^{\infty}$ with
$(f_{\theta}/f_{\theta_{0}}(k))_{k=1}^{\infty}$.

Fix any $\theta\in\Theta^{\prime}$ and assume that $\theta$ renders
$f_{\theta}/f_{\theta_{0}}$ decreasing. By Lemma 4.4 applied with weights
$w_{k}=f_{\theta}/f_{\theta_{0}}(k)$, there exists $K$ such that
$p_{>k}^{\scriptscriptstyle(0)}\leq q^{\scriptscriptstyle(0,\theta)}_{k}$ for
$k\leq K$ and
$p_{>k}^{\scriptscriptstyle(0)}>q^{\scriptscriptstyle(0,\theta)}_{k}$ for
$k>K$. Then Lemma 4.3 with $v_{k}=\dot{f}_{\theta}/f_{\theta}(k)$ (understood
component-wise) and the probability distributions
$p_{>k}^{\scriptscriptstyle(0)}$ and $q_{k}^{\scriptscriptstyle(0,\theta)}$,
shows that $\dot{\iota}(f_{\theta})<0$. If $f_{\theta}/f_{\theta_{0}}(k)$ is
increasing, then the same argument applies, but we find that
$\dot{\iota}(f_{\theta})>0$.

For reference, we state the martingale central limit theorem—a version of
Theorem 3.2 of Hall and Heyde (2014).

######  Proposition 9.2.

Suppose that $X_{t}$ is a martingale difference series relative to the
filtration $\mathcal{F}_{t}$. If as $n\rightarrow\infty$,
$n^{-1}\sum_{t=1}^{n}\mathbb{E}[X_{t}^{2}|\mathcal{F}_{t-1}]\xrightarrow{P}v$
for a positive constant $v$ and
$n^{-1}\sum_{t=1}^{n}\mathbb{E}[X^{2}_{t}{1}_{\\{|X_{t}|>\varepsilon\sqrt{n}\\}}|\mathcal{F}_{t-1}]\xrightarrow{P}0$
for every $\varepsilon>0$, then $\sqrt{n}\bar{X}_{n}\rightsquigarrow N(0,v)$.

###### Proof of Theorem 5.1.

By Theorem 4.1 the maximum likelihood estimator $\hat{\theta}_{n}$ is
consistent if ${\theta}_{0}$ is identifiable. The $\hat{\theta}_{n}$ will
eventually be interior to the parameter set if ${\theta}_{0}$ is interior, and
hence satisfy the system of likelihood equations
$\dot{\iota}_{n}(f_{\hat{\theta}_{n}})=0$. Thus, the second assertion of the
theorem follows from the first.

By a Taylor expansion of the $i$th of the likelihood equations, it can be
expanded as
$0=\dot{\iota}_{n}(f_{\theta_{0}})_{i}+\ddot{\iota}_{n}(f_{\theta^{\prime}_{n,i}})_{i}(\hat{\theta}_{n}-\theta_{0})$,
for $\theta^{\prime}_{n,i}$ on the line segment between $\hat{\theta}_{n}$ and
$\theta_{0}$ and $\ddot{\iota}_{n}(f_{\theta})_{i}$ the $i$th row of the
second derivative matrix $\ddot{\iota}_{n}(f_{\theta})$. Thus,
$\ddot{\iota}_{n}(f_{\theta^{\prime}_{n}})(\hat{\theta}_{n}-\theta_{0})=-\dot{\iota}_{n}(f_{\theta_{0}})$,
where $\ddot{\iota}_{n}(f_{\theta^{\prime}_{n}})$ is understood to be the
$(d\times d)$ matrix with $i$th row
$\ddot{\iota}_{n}(f_{\theta^{\prime}_{n,i}})_{i}$, even though the vector
$\theta^{\prime}_{n,i}$ may be different for different $i$. The proof can be
concluded by showing that
$\sqrt{n}\,\dot{\iota}_{n}(f_{\theta_{0}})\rightsquigarrow N(0,V_{0})$ and
that $\ddot{\iota}_{n}(f_{\theta^{\prime}_{n}})\rightarrow-V_{0}$ in
probability.

As noted in Section 2, the sequence $n\dot{\iota}_{n}(f_{\theta_{0}})$ is a
(vector-valued) martingale. The asymptotic normality can be obtained from the
martingale central limit theorem (see Proposition 9.2). Because
$\mathbb{P}_{\theta}(D_{t}=k\mathchar 25194\relax{\cal
F}_{t-1})=f_{{\theta}}(k)N_{k}(t-1)/S_{f_{\theta}}(t-1)$, the martingale
differences
$(\dot{f}_{{\theta}_{0}}/f_{{\theta}_{0}})(D_{t})-\mathbb{E}_{{\theta}_{0}}\bigl{(}(\dot{f}_{{\theta}_{0}}/f_{{\theta}_{0}})(D_{t})\mathchar
25194\relax{\cal F}_{t-1}\bigr{)}$ possess conditional covariances

$\displaystyle\Sigma_{t}:=$
$\displaystyle\sum_{k=1}^{\infty}\Bigl{(}\frac{\dot{f}_{\theta_{0}}}{f_{\theta_{0}}}\Bigr{)}\Bigl{(}\frac{\dot{f}_{\theta_{0}}}{f_{\theta_{0}}}\Bigr{)}^{T}(k)\frac{f_{\theta_{0}}(k)P_{k}(t-1)}{S_{f_{{\theta}_{0}}}(t-1)/(t-1)}$
$\displaystyle\qquad-\Bigl{(}\sum_{k=1}^{\infty}\frac{\dot{f}_{\theta_{0}}}{f_{\theta_{0}}}(k)\frac{f_{\theta_{0}}(k)P_{k}(t-1)}{S_{f_{{\theta}_{0}}}(t-1)/(t-1)}\Bigr{)}\Bigl{(}\sum_{k=1}^{\infty}\frac{\dot{f}_{\theta_{0}}}{f_{\theta_{0}}}(k)\frac{f_{\theta_{0}}(k)P_{k}(t-1)}{S_{f_{{\theta}_{0}}}(t-1)/(t-1)}\Bigr{)}^{T}.$

As seen in the proof of Theorem 4.1, the sequence $S_{f_{{\theta}_{0}}}(t)/t$
tends almost surely to
$\sum_{j}f_{{\theta}_{0}}(j)p_{j}^{\scriptscriptstyle(0)}$, as
$t\rightarrow\infty$. Corollary 3.2 applied with $h$ equal to the entries of
the matrix $\dot{f}_{{\theta}_{0}}\dot{f}_{{\theta}_{0}}^{T}/f_{{\theta}_{0}}$
or the vector $\dot{f}_{{\theta}_{0}}$ shows that the preceding display tends
almost surely to

$\sum_{k=1}^{\infty}\frac{\dot{f}_{\theta_{0}}\dot{f}_{{\theta}_{0}}^{T}}{f_{\theta_{0}}}(k)\frac{p_{k}^{\scriptscriptstyle(0)}}{\sum_{j}f_{{\theta}_{0}}(j)p_{j}^{\scriptscriptstyle(0)}}-\Bigl{(}\sum_{k=1}^{\infty}\dot{f}_{\theta_{0}}(k)\frac{p_{k}^{\scriptscriptstyle(0)}}{\sum_{j}f_{{\theta}_{0}}(j)p_{j}^{\scriptscriptstyle(0)}}\Bigr{)}\Bigl{(}\sum_{k=1}^{\infty}\dot{f}_{\theta_{0}}(k)\frac{p_{k}^{\scriptscriptstyle(0)}}{\sum_{j}f_{{\theta}_{0}}(j)p_{j}^{\scriptscriptstyle(0)}}\Bigr{)}^{T}.$

In view of equation (3.8), this is equal to the matrix $V_{0}$. The averages
$n^{-1}\sum_{t=2}^{n}\Sigma_{t}$ of the conditional covariances tend to the
same limit, by Lemma 9.1.

Because $D_{t}\leq t$, we bound, using (5.2),

$\Bigl{\|}\frac{\dot{f}_{{\theta}_{0}}}{f_{{\theta}_{0}}}(D_{t})\Bigr{\|}\leq
C\log^{\gamma}t\leq C\log^{\gamma}n,\qquad t\leq n.$

As this is smaller than $\epsilon\sqrt{n}$, eventually for every $\epsilon>0$,
the conditional Lindeberg condition is trivially satisfied. We conclude that
$\sqrt{n}\,\dot{\iota}_{n}(f_{\theta_{0}})\rightsquigarrow N(0,V_{0})$, by the
martingale central limit theorem, for instance, Proposition 9.2.

The Hessian matrix $\ddot{\iota}_{n}(\theta)$ takes the form

$\displaystyle\ddot{\iota}_{n}(f_{\theta})$
$\displaystyle=\sum_{k=1}^{\infty}\Bigl{(}\frac{\ddot{f}_{\theta}}{f_{\theta}}-\frac{\dot{f}_{\theta}\dot{f}_{\theta}^{T}}{f_{\theta}^{2}}\Bigr{)}(k)P_{>k}(n)-\frac{1}{n}\sum_{t=2}^{n}\Bigl{(}\frac{S_{\ddot{f}_{\theta}}}{S_{f_{\theta}}}-\frac{S_{\dot{f}_{\theta}}S_{\dot{f}_{\theta}^{T}}}{S_{\dot{f}_{\theta}}^{2}}\Bigr{)}(t-1)$
$\displaystyle=\sum_{j=1}^{\infty}\sum_{k=1}^{j-1}\Bigl{(}\frac{\ddot{f}_{\theta}}{f_{\theta}}-\frac{\dot{f}_{\theta}\dot{f}_{\theta}^{T}}{f_{\theta}^{2}}\Bigr{)}(k)P_{j}(n)-\frac{1}{n}\sum_{t=2}^{n}\Bigl{(}\frac{S_{\ddot{f}_{\theta}}}{S_{f_{\theta}}}-\frac{S_{\dot{f}_{\theta}}S_{\dot{f}_{\theta}^{T}}}{S_{\dot{f}_{\theta}}^{2}}\Bigr{)}(t-1).$
(9.3)

Using Corollary 3.2, the first term in the second line can be shown to
converge as $n\rightarrow\infty$ to the expression obtained by replacing
$P_{j}(n)$ by $p_{j}^{\scriptscriptstyle(0)}$, or equivalently replacing
$P_{>k}(n)$ by $p_{>k}^{\scriptscriptstyle(0)}$ in the first line. For the
second term we first note that $S_{f}(t)/t$ tends almost surely to
$\sum_{j}f(j)p_{j}^{\scriptscriptstyle(0)}$, as $t\rightarrow\infty$, for $f$
equal to $\ddot{f}_{{\theta}}$, $\dot{f}_{\theta}$ or $f_{\theta}$. By the
continuous mapping theorem the terms of the sum converge to the corresponding
limit. The second term then converges almost surely to the same limit, in view
of Lemma 9.1, still uniformly in ${\theta}$. By arguments similar to those in
the proof of Theorem 4.1, the convergences of both terms can be seen to be
uniform in ${\theta}$.

Finally, the continuity of the limit and consistency of
$\hat{\theta}_{n}^{\prime}$ for ${\theta}_{0}$ give that the
$\ddot{\iota}_{n}(f_{\theta_{n}^{\prime}})$ tends to the limit evaluated at
${\theta}_{0}$. This can be seen to be equal to the matrix $-V_{0}$ with the
help of (3.8), where the two terms involving $\ddot{f}_{\theta}$ cancel each
other.

###### Proof of Theorem 5.3.

Because the true parameter is on the boundary of the parameter set, the
maximum likelihood estimator may not solve the likelihood equations. Instead,
we use its characterization as the maximizer of the log likelihood. The log
likelihood evaluated at the parameter ${\theta}_{0}+h/\sqrt{n}$ satisfies

$\ell_{n}\bigl{(}f_{{\theta}_{0}+h/\sqrt{n}}\bigr{)}-\ell_{n}(f_{{\theta}_{0}})=h^{T}\sqrt{n}\,\dot{\iota}_{n}(f_{{\theta}_{0}})+{\textstyle{\frac{1}{2}}}h^{T}\ddot{\iota}_{n}(f_{{\theta}_{n}^{\prime}(h)})h,$

where ${\theta}_{n}^{\prime}(h)$ is on the line segment between ${\theta}_{0}$
and ${\theta}_{0}+h/\sqrt{n}$. The rescaled maximum likelihood estimator
$\hat{h}_{n}=\sqrt{n}(\hat{\theta}_{n}-{\theta}_{0})$ maximizes this process
over the set $H_{n}$ of all $h=\sqrt{n}({\theta}-{\theta}_{0})$ such that
${\theta}=({\alpha},{\beta})$ belongs to the parameter set
$\\{(a,{\beta}):{\alpha}\in(-1+{\varepsilon},{\varepsilon}^{-1}),{\beta}\leq
1\\}$. Since $\hat{\theta}_{n}$ is consistent for ${\theta}_{0}$, by Theorem
4.1, the set $H_{n}$ can be further reduced to a set such that
$\|h\|<\sqrt{n}\,{\delta}_{n}$, for some ${\delta}_{n}\rightarrow 0$. By the
arguments in the proof of Theorem 5.1, we have
$\sup_{\theta}|\ddot{\iota}_{n}(f_{\theta})-\ddot{\iota}(f_{\theta})|\rightarrow
0$, in probability. Combined with the continuity of
${\theta}\mapsto\ddot{\iota}(f_{\theta})$, we see that $\sup_{h\in
H_{n}}|\ddot{\iota}_{n}(f_{{\theta}_{n}^{\prime}(h)})-\ddot{\iota}(f_{{\theta}_{0}})|\rightarrow
0$, almost surely. From the nonsingularity and negative definiteness of
$-V_{0}=\ddot{\iota}(f_{{\theta}_{0}})$, we then see that
$h^{T}\ddot{\iota}_{n}(f_{{\theta}_{n}^{\prime}(h)})h<-c\|h\|^{2}$, for every
$h\in H_{n}$ and some $c>0$, with probability tending to one. Using that $0\in
H_{n}$ and
$|h^{T}\sqrt{n}\,\dot{\iota}_{n}(f_{{\theta}_{0}})|\leq\|h\|O_{P}(1)$, we
conclude that $\hat{h}_{n}=O_{P}(1)$. Next the argmax continuous mapping
theorem (e.g. Corollary 5.58 and Lemma 7.13 in van der Vaart (2000)) shows
that

$\hat{h}_{n}=\mathop{\rm argmax}_{h\in
H_{n}}\bigl{(}h^{T}\sqrt{n}\,\dot{\iota}_{n}(f_{{\theta}_{0}})+{\textstyle{\frac{1}{2}}}h^{T}\ddot{\iota}_{n}(f_{{\theta}_{n}^{\prime}(h)})h\bigr{)}\rightsquigarrow\mathop{\rm
argmax}_{h\in H}(h^{T}Z_{0}-{\textstyle{\frac{1}{2}}}h^{T}V_{0}h),$

for $H=\\{(h_{1},h_{2}):h_{2}\leq 0\\}$ the limit of the sequence of sets
$H_{n}$ and $Z_{0}\sim N(0,V_{0})$ the limit in distribution of the sequence
$\sqrt{n}\,\dot{\iota}_{n}(f_{{\theta}_{0}})$. The right side has the claimed
distribution, by Lemma 9.3 (where we set $V=V_{0}$).

######  Lemma 9.3.

If $\hat{h}=\mathop{\rm argmax}_{h:a^{T}h\leq 0}(2h^{T}Z-h^{T}V_{0}h)$ for
$a\in\mathbb{R}^{d}$ and a random variable $Z\sim N_{d}(0,V)$ for positive
semi-definite matrices $V$ and $V_{0}$, then $a^{T}\hat{h}\sim W{1}_{W\leq 0}$
for $W\sim N(0,a^{T}V_{0}^{-1}VV_{0}^{-1}a)$.

###### Proof of Lemma 9.3.

Since
$2h^{T}Z-h^{T}V_{0}h=-\|V_{0}^{1/2}h-V_{0}^{-1/2}Z\|^{2}+Z^{T}V_{0}^{-1}Z$,
the variable $\hat{h}$ can be seen to be equal to
$V_{0}^{-1/2}\Pi_{G}(V_{0}^{-1/2}Z)$, for $\Pi_{G}$ the projection onto the
half space $G:=V_{0}^{1/2}\\{h:a^{T}h\leq 0\\}=\\{g:b^{T}g\leq 0\\}$, for
$b=V_{0}^{-1/2}a$. Therefore,
$a^{T}\hat{h}=b^{T}\Pi_{G}(\tilde{W})={1}_{b^{T}\tilde{W}\leq
0}b^{T}\tilde{W}$, for $\tilde{W}:=V_{0}^{-1/2}Z\sim
N_{d}(0,V_{0}^{-1/2}VV_{0}^{-1/2})$. The proof is complete by setting
$W:=b^{T}\tilde{W}=a^{T}V_{0}^{-1}Z$.

###### Proof of Theorem 6.2.

Since
$\sum_{k=1}^{\infty}a_{k}P_{>k}(n)=\sum_{j=1}^{\infty}\sum_{k=1}^{j-1}a_{k}P_{j}(n)$,
we can write $\dot{\tilde{\iota}}_{n}(f_{\theta_{0}})=\phi\bigl{(}S_{n})$, for
$S_{n}=\sum_{k=1}^{\infty}F_{{\theta}_{0}}(k)P_{k}(n)$, and
$\phi:\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{R}\to\mathbb{R}^{d}$
given by $\phi(s_{1},s_{2},s_{3})=s_{1}-\frac{1}{s_{3}}s_{2}$. Similarly, it
follows from (3.8) that $\phi(S^{\scriptscriptstyle(0)})=0$, for
$S^{\scriptscriptstyle(0)}=\sum_{k=1}^{\infty}F_{{\theta}_{0}}(k)p_{k}^{\scriptscriptstyle(0)}$.
The assumption on $S_{n}$ and the delta-method then give that the sequence
$\sqrt{n}\dot{\tilde{\iota}}_{n}(f_{\theta_{0}})-\sqrt{n}\bigl{(}\phi(S_{n})-\phi(S^{\scriptscriptstyle(0)})\bigr{)}$
tends in $\mathbb{R}^{d}$ in distribution to a normal distribution with mean
zero and covariance matrix
$V=(I,-I/S^{\scriptscriptstyle(0)}_{3},S^{\scriptscriptstyle(0)}_{2}/(S^{\scriptscriptstyle(0)}_{3})^{2})W_{0}(I,-I/S^{\scriptscriptstyle(0)}_{3},S^{\scriptscriptstyle(0)}_{2}/(S^{\scriptscriptstyle(0)}_{3})^{2})^{T}$.

The difference between the pseudo score function
$\dot{\tilde{\iota}}(f_{\theta})$ and $\dot{\iota}(f_{\theta})$ is only in
their second terms, which are $S_{\dot{f}_{\theta}}/S_{f_{\theta}}(n)$ for the
first and the Césaro averages of these for the second. The Hessian matrices
$\ddot{\tilde{\iota}}(f_{\theta})$ and $\ddot{\iota}(f_{\theta})$ also only
differ by the derivatives of the second terms. In the proof of Theorem 5.1 it
was noted that the latter are again Césaro averages (see (9.3)) and they were
shown to converge by showing convergence of the individual terms. Thus, this
proof applies also in the present situation, and the full proof can be
finished as the proof of Theorem 5.1.

###### Proof of Corollary 6.3.

In view of Theorem 6.2 it suffices to verify the asymptotic normality of the
sequence
$\sum_{k=1}^{\infty}F_{{\theta}_{0}}(k)\sqrt{n}(P_{k}(n)-p_{k}^{\scriptscriptstyle(0)})$,
for $F_{{\theta}_{0}}$ as given. Because $f_{{\theta}}(k)$ is constant in
$k\geq K$, so is $F_{{\theta}_{0}}(k)$ and hence
$\sum_{k}F_{{\theta}_{0}}(k)P_{k}(n)=\sum_{k\leq
K}F_{{\theta}_{0}}(k)P_{k}(n)+F_{{\theta}_{0}}(K)P_{>K}(n)$, and the same for
$p_{k}^{\scriptscriptstyle(0)}$ instead of $P_{k}(n)$. The desired convergence
therefore follows from Proposition B.4 and the continuous mapping theorem.

##  Appendix A Branching Processes and Rooted Ordered Trees

A rooted ordered tree is a tree in which one node is designated as the root
and the other nodes can be oriented in parent-child relations in reference to
their distance to this root node. In a dynamic setup the root is the initial
ancestor who is responsible for giving births directly or indirectly to every
other node. In the Ulam–Harris labelling notation for branching processes, the
root is denoted by $\emptyset$ and every other node has the form
$(i_{1},\dots,i_{l})$, for positive natural numbers $i_{j}\in\mathbb{N}_{+}$
with $l\in\mathbb{N}_{+}$. The node $x=(i)$ is the $i$-th child of the root,
and more generally the node $x=(i_{1},\dots,i_{k})$ is the $i_{k}$-th child of
$\tilde{x}=(i_{1},\dots,i_{k-1})$. By induction, the set of all possible
individuals is

$\mathscr{I}=\\{\emptyset\\}\cup\left(\bigcup_{k=1}^{\infty}\mathbb{N}_{+}^{k}\right).$

The root is the zero-th generation and the $k$-th generation consists of all
$x\in\mathbb{N}_{+}^{k}$. For $x=(x_{1},\ldots,x_{k})$ and
$y=(y_{1},\dots,y_{l})$ the notation $xy$ is shorthand for the concatenation
$(x_{1},\dots,x_{k},y_{1},\dots,y_{l})$, and, in particular,
$xl=(x_{1},\dots,x_{k},l)$. The labelling of the nodes contains all parental
information and a rooted ordered tree is defined to be a subset
$G\subset\mathscr{I}$ such that if $x=(x_{1},\ldots,x_{k})\in G$, then both
$(x_{1},\ldots,x_{k-1})\in G$ and $(x_{1},\ldots,x_{k}-1)\in G$ in case
$x_{k}\geq 2$. A corresponding graphical representation is obtained by drawing
a node for every $x\in G$ and connecting two nodes by an edge if they are in a
parent-child relationship. The degree of node $x\in G$ is

$\deg(x,G)=|\\{l\in\mathbb{N}_{+}\mid xl\in G\\}|+1,$ (A.1)

where the extra one is for the parent of the node.

To set up a stochastic branching process, each individual $x\in\mathscr{I}$ is
associated with a stochastic variable $\lambda_{x}$ and two stochastic
processes $\xi_{x}$ and $\varphi_{x}$, called the life length, the
reproduction process and the characteristic of $x$. The triples
$(\lambda_{x},\xi_{x},\varphi_{x})$ are taken IID across the nodes
$x\in\mathscr{I}$. Formally they may be defined as copies on a product
probability space

$(\Omega,\mathcal{B},P)=\prod_{x\in\mathscr{I}}(\Omega_{x},\mathcal{B}_{x},P_{x}),$

where each $(\Omega_{x},\mathcal{B}_{x},P_{x})$ is a copy of a probability
space $(\Omega_{0},\mathcal{B}_{0},P_{0})$. For a given measurable map
$(\lambda,\xi,\varphi)$ defined on $(\Omega_{0},\mathcal{B}_{0},P_{0})$, we
then define
$(\lambda_{x},\xi_{x},\varphi_{x})(\omega)=(\lambda,\xi,\varphi)(\omega_{x})$
if $\omega=(\omega_{x})_{x\in\mathscr{I}}\in\Omega$. The life length $\lambda$
is a nonnegative random variable, which in our case we take identically
$\infty$ (no node will die). The reproduction process $\xi=(\xi(t):t\geq 0)$
will be a counting process starting with $\xi(0)=0$ and increasing by steps of
size $1$ at random times. We identify $\xi$ with a random
$\mathbb{N}_{+}$-valued measure through $\xi([0,t])=\xi(t)$ and denote by
$\mu:t\mapsto\mathbb{E}[\xi(t)]$ its mean (or intensity) measure, which is
often called the reproduction function in this context. The characteristic
will also be a stochastic process $\varphi=(\varphi(t):t\geq 0)$, where we may
set $\varphi(t)=0$ for $t<0$.

The random point process $\xi_{x}$ models the birth times of the children of
individual $x$ relative to the birth time $\sigma_{x}$ of $x$. The latter
birth times are formally defined recursively by setting the birth time of the
root $\emptyset$ at $t=0$ (hence $\sigma_{\emptyset}=0$), and next the birth
time $\sigma_{y}$ of $y$ in calendar time by

$\sigma_{y}=\sigma_{x}+\inf\\{t\geq 0\colon\xi_{x}(t)\geq l\\},\qquad\text{ if
}y=xl.$

The calendar time is the evolution time of the branching process, as opposed
to the local time scales of the processes $\xi_{x}$ and $\varphi_{x}$, of
which the local zero time is interpreted to be $\sigma_{x}$ in calendar time.
The variable $\varphi_{x}(t)$ is interpreted as the characteristic of
individual $x$ when $x$ has age $t$. Calendar time is also different from the
discrete time steps used to describe the evolution of a PA network.

For a given characteristic we define the process

$Z_{t}^{\varphi}=\sum_{x\in\mathscr{I}\colon\sigma_{x}\leq
t}\varphi_{x}(t-\sigma_{x}).$

The variable $t-\sigma_{x}$ is the time since birth of individual $x$ and
hence $\varphi_{x}(t-\sigma_{x})$ can be interpreted as the characteristic of
individual $x$ at calendar time $t$. The variable $Z_{t}^{\varphi}$ is the sum
of all such characteristics over the individuals that are alive at time $t$.

A branching process is supercritical and Malthusian if its reproduction
function $\mu$ does not concentrate on any lattice $\\{0,h,2h,\dots\\}$, for
some $h>0$, and there exists a number ${\lambda^{*}}>0$ such that

$\int_{0}^{\infty}{\mathrm{e}}^{-{\lambda^{*}}t}\mu(dt)=1.$ (A.2)

We shall also assume the integrability assumption

$\int_{0}^{\infty}t^{2}\mathrm{e}^{-{\lambda^{*}}t}\,\mu(dt)<\infty.$ (A.3)

Existence of a solution to equation (A.2) is called the Malthusian assumption,
and ${\lambda^{*}}$ is called the Malthusian parameter.

The following proposition can be obtained by combining Theorem 3.1, Corollary
3.4 and Theorem 6.3 of Nerman (1981). Define
${}_{\lambda}\xi(t)=\int_{0}^{t}e^{-\lambda u}\,\xi(du)$.

######  Proposition A.1.

Assume that the reproduction function $\mu(t)=\mathbb{E}[\xi(t)]$ satisfies
conditions (A.2) and (A.3) and does not concentrate on any lattice. Assume
that $t\mapsto\mathbb{E}[\varphi(t)]$ is continuous almost everywhere with
respect to the Lebesgue measure and the following conditions hold:

$\displaystyle\sum_{k=0}^{\infty}\sup_{k\leq t\leq
k+1}\big{(}e^{-{\lambda^{*}}t}\mathbb{E}[\varphi(t)]\big{)}<\infty,$ (A.4)
$\displaystyle\mathbb{E}[\sup_{s\leq t}\varphi(s)]<\infty\quad\text{for all
}t<\infty.$ (A.5)

Then there exists a random variable $Y_{\infty}$ depending only on the
reproduction process $\xi(t)$ such that, as $t\rightarrow\infty$,

$e^{-{\lambda^{*}}t}Z_{t}^{\varphi}\xrightarrow{P}Y_{\infty}m_{\infty}^{\varphi},$
(A.6)

where $m_{\infty}^{\varphi}$ is defined as

$m_{\infty}^{\varphi}=\frac{\int_{0}^{\infty}e^{-{\lambda^{*}}t}\mathbb{E}[\varphi(t)]\,dt}{\int_{0}^{\infty}te^{-{\lambda^{*}}t}\,d\mu(t)}.$

The convergence in (A.6) also holds in the $L_{1}$ sense if

$\mathbb{E}[_{\lambda^{*}}\xi(\infty)\log^{+}{}_{\lambda^{*}}\xi(\infty)]<\infty.$
(A.7)

Suppose that the reproduction process $\xi$ satisfies (A.7), and both
$\varphi_{1}$ and $\varphi_{2}$ satisfy the conditions (A.4) and (A.5). Define
$T_{t}$ as total number of births up to and including time $t$. Then, on the
event $\\{T_{t}\rightarrow\infty\\}$, as $t\rightarrow\infty$.

$\frac{Z_{t}^{\varphi_{1}}}{Z_{t}^{\varphi_{2}}}\xrightarrow{P}\frac{m_{\infty}^{\varphi_{1}}}{m_{\infty}^{\varphi_{2}}}=\frac{\int_{0}^{\infty}e^{-{\lambda^{*}}t}\mathbb{E}[\varphi_{1}(t)]\,dt}{\int_{0}^{\infty}e^{-{\lambda^{*}}t}\mathbb{E}[\varphi_{2}(t)]\,dt}.$
(A.8)

If $\varphi_{1}$ and $\varphi_{2}$ have càdlàg paths and there exists a
$\lambda<{\lambda^{*}}$ such that

$\displaystyle\mathbb{E}[_{\lambda}\xi(\infty)]$ $\displaystyle<\infty,$ (A.9)
$\displaystyle\mathbb{E}\bigl{[}\sup_{t}e^{-\lambda t}\varphi_{i}(t)\bigr{]}$
$\displaystyle<\infty,\qquad i=1,2,$ (A.10)

then, on $\\{T_{t}\rightarrow\infty\\}$, the convergence in (A.8) is also in
the almost sure sense.

##  Appendix B Asymptotic normality of Empirical Degrees

In this section we derive the asymptotic normality of the empirical degrees in
some cases of the preferential attachment model, using an urn process studied
by Janson (2004).

The urn process consists of vectors $X_{n}=(X_{n,1},\dots,X_{n,q})^{T}$ in
$[0,\infty)^{q}$, of which the $i$-th coordinate represents the quantity in
the $i$-th urn at time $n$. Given are, for each $i\in[q]$, an ‘activity’
$a_{i}\geq 0$ and a vector
$\xi_{i}=(\xi_{i,1},\dots,\xi_{i,q})^{T}\in\mathbb{R}^{q}$ with $\xi_{i,j}\geq
0$ for $j\neq i$ and $\xi_{ii}\geq-1$. The process $(X_{n})_{n=0}^{\infty}$
evolves as a Markov process, with transitions determined by: given $X_{n-1}$,

  1. 1.

pick an urn $i\in[q]$ with probability
$a_{i}X_{n-1,i}/\sum_{j=1}^{q}a_{j}X_{n-1,j}$;

  2. 2.

set $X_{n}:=X_{n-1}+\xi_{i}$.

It is assumed that $\sum_{j=1}^{q}\xi_{i,j}\geq 0$, for every $i\in[q]$, with
strict inequality for some $i$, so that the total content
$\sum_{i=1}^{q}X_{n,i}$ of the urns is nondecreasing. (Actually, Janson (2004)
allows the vectors $\xi_{i}$ to be random, but deterministic vectors suffice
in our situation, and allow a simpler statement of the main result.)

Define the ‘transfer’ matrix $A\in\mathbb{R}^{q\times q}$ by

$A_{ij}=a_{j}\xi_{j,i}.$ (B.1)

We assume that $A$ is irreducible. By an application of the Perron–Frobenius
theorem (to the nonnegative matrix $A+{\alpha}I_{q}$, for sufficiently large
${\alpha}$), it can be seen that the eigenvalue of $A$ with the largest real
value is real, and the eigenvalues can be ranked by their real parts as
$\lambda_{1}>\mathop{\rm Re}\nolimits\lambda_{2}\geq\mathop{\rm
Re}\nolimits\lambda_{3}\geq\cdots$. Furthermore, ${\lambda}_{1}>0$, has
multiplicity one, and the associated eigenvector $v_{1}$ has all positive
coordinates. We normalize $v_{1}$ such that $a^{T}v_{1}=1$, where
$a=(a_{1},\ldots,a_{q})^{T}$.

In this setup, conditions (A1)-(A6) in Janson (2004) are satisfied (see his
Lemma 2.1). A crucial further assumption in the following proposition, which
restates Theorems 3.21-3.22 of Janson (2004), is that $\mathop{\rm
Re}\nolimits\lambda_{2}<{\lambda}_{1}/2$.

Define the following quantities:

$\displaystyle B:=\sum_{i=1}^{q}v_{1i}a_{i}\xi_{i}\xi_{i}^{T},$ (B.2)
$\displaystyle\varphi(s,A):=\sum_{n=1}^{\infty}\frac{s^{n}}{n!}A^{n-1}=\int_{0}^{s}e^{tA}\,dt,$
(B.3) $\displaystyle\psi(s,A):=e^{sA}-\lambda_{1}v_{1}a^{T}\varphi(s,A).$
(B.4)

######  Proposition B.1.

Under the preceding conditions, as $n\rightarrow\infty$,

$n^{-1}X_{n}\xrightarrow{\text{a.s.}}\lambda_{1}v_{1}.$

If, moreover, $\mathop{\rm Re}\nolimits\lambda_{2}<\lambda_{1}/2$, then, as
$n\rightarrow\infty$,

$n^{1/2}(n^{-1}X_{n}-\lambda_{1}v_{1})\xrightarrow{\text{d}}N(0,\Sigma),$

where the covariance matrix $\Sigma$ is defined as (with quantities defined in
(B.2)–(B.4))

$\Sigma=\int_{0}^{\infty}\psi(s,A)B\psi(s,A)^{T}e^{-\lambda_{1}s}\lambda_{1}ds-\lambda_{1}^{2}v_{1}v_{1}^{T}.$
(B.5)

The preferential attachment model would be naturally described using an
infinite number of urns, with the content $X_{n,i}$ of the $i$-th urn
corresponding to the number of nodes of degree $i$ at time $n$. The activities
$a_{i}$ can then be set equal to the preferences $f(i)$ and the vectors
$\xi_{i}$ defined by their coordinates

$\xi_{i,j}=-{1}_{\\{j=i\\}}+{1}_{\\{j=i+1\\}}+{1}_{\\{j=1\\}}.$ (B.6)

The last term ${1}_{\\{j=1\\}}$ corresponds to the new node, which has degree
1, and is counted in the first urn, while the first two terms on the right
describe the decrease and increase by 1 of the numbers of nodes of degrees $j$
and $j+1$, respectively, if the new node is attached to an existing node of
degree $j$. There would then be infinitely many vectors ($i\in\mathbb{N}_{+}$)
each with infinitely many coordinates ($j\in\mathbb{N}_{+}$), but,
unfortunately, Proposition B.5 allows a fixed, finite number of urns only. In
the following, we consider two examples of PA models in which the infinite
process can be reduced to a finite number: the case of an affine PA function,
and the case that the PA function is constant from a fixed degree onward.

In the affine case with PA function $f(k)=k+{\alpha}$, we study the degree
distribution up to some given degree ${\kappa}$ by gathering all nodes of
degree strictly bigger than ${\kappa}$ in a single urn, labelled
${\kappa}+1=:q$. To accommodate that the latter nodes have different
preferences, we define the vectors $\xi_{i}\in\mathbb{R}^{{\kappa}+1}$ for
$i=1,\ldots,{\kappa}-1$ as in (B.6) with coordinates restricted to
$j\in[{\kappa}+1]$, but redefine $\xi_{\kappa}$ and $\xi_{{\kappa}+1}$ by

$\displaystyle\xi_{\kappa,j}$
$\displaystyle=-{1}_{\\{j=\kappa\\}}+{1}_{\\{j=\kappa+1\\}}(\kappa+1+{\alpha})+{1}_{\\{j=1\\}},$
(B.7) $\displaystyle\xi_{\kappa+1,j}$
$\displaystyle={1}_{\\{j=\kappa+1\\}}+{1}_{\\{j=1\\}}.$ (B.8)

We combine this with the vector of activities
$a_{\kappa}=(1+{\alpha},2+{\alpha},\dots,\kappa+{\alpha},1)^{T}$. Thus, urns
$1,\ldots,{\kappa}$ have activities equal to the preferential attachment
function, but urn ${\kappa}+1$ has activity 1. With these definitions, for
$i=1,\ldots,{\kappa}$ the variable $X_{n,i}$ corresponds to the numbers of
nodes of degree $i$, but $X_{n,{\kappa}+1}$ is set to correspond to the total
preference of all nodes of degree bigger than ${\kappa}$. Indeed, a choice of
an urn $i=1,\ldots,{\kappa}-1$ follows the scheme described before, with the
transition given by (B.6). Second, a choice of urn ${\kappa}$ corresponds to
choosing a node of degree ${\kappa}$; by (B.7) the count of urn ${\kappa}$ is
then decreased by one, and ${\kappa}+1+{\alpha}$ balls are added to urn
${\kappa}+1$, each weighted by activity $a_{{\kappa}+1}=1$, thus giving the
correct increase of total preference of the nodes of degree bigger than
${\kappa}$. Third, a choice of urn ${\kappa}+1$ corresponds to choosing a node
of some degree bigger than ${\kappa}$; this node is replaced by a node of
degree one bigger, resulting in an increase by 1 of the total preference of
the nodes of degree bigger than ${\kappa}$. Thus, (B.8) correctly changes the
total preferences of the nodes of degree bigger than ${\kappa}$. In both (B.7)
and (B.8) the term on the far right corresponds to the new node of degree 1,
counted in $X_{n,1}$.

Asymptotic normality of the empirical degrees $P_{k}(n)$ in the affine case
was first proved in Móri (2002) and Resnick and Samorodnitsky (2016). We
deduce it here by a simple argument based on the preceding proposition,

######  Proposition B.2.

In the PA model with PA function $f(k)=k+\alpha$, the centered and rescaled
empirical degree distribution
$\bigl{(}\sqrt{n}\bigl{(}P_{k}(n)-p_{k}\bigr{)},k=1,2,\dots\bigr{)}$ converges
in distribution in $\mathbb{R}^{\mathbb{N}_{+}}$ to a centered Gaussian
process.

###### A novel proof.

We fix arbitrary ${\kappa}\in\mathbb{N}$, and define activities $a_{i}$ and
update vectors ${\xi}_{i}$, for $i\in[{\kappa}+1]$, as indicated in
(B.6)–(B.8), and initial vector $X_{0}=(1,0,\dots,0)^{T}$. For every $n$, the
vector $(X_{n,1},\ldots,X_{n,{\kappa}},X_{n,{\kappa}+1})$ is then identically
distributed to the vector
$\bigl{(}N_{1}(n),\ldots,N_{\kappa}(n),\sum_{j>{\kappa}}N_{j}(n)(j+{\alpha})\bigr{)}$.
Since (see, e.g., Billingsley (2013, Example 2.4) or van der Vaart and Wellner
(1996, Theorem 1.6.1)) weak convergence in $\mathbb{R}^{\mathbb{N}_{+}}$ is
the same as convergence of all finite marginals, the first assertion is proved
if we can prove convergence of these vectors for every fixed ${\kappa}$. For
this we apply Proposition B.2.

The transfer matrix $A=A_{\kappa}$ is given by

$A_{\kappa}=\begin{pmatrix}0&\ 2+{\alpha}&\
3+{\alpha}&\cdots&\kappa-1+{\alpha}&\kappa+{\alpha}&1\\\ 1+{\alpha}\ \
&-2-{\alpha}&0&\cdots&0&0&0\\\ 0&\ 2+{\alpha}&-3-{\alpha}\ \ &\
\cdots&0&0&0\\\ 0&0&\ 3+{\alpha}&\cdots&0&0&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
0&0&0&\cdots&-{\kappa}+1-{\alpha}&0&0\\\ 0&0&0&\cdots&\
{\kappa}-1+{\alpha}&-\kappa-{\alpha}&0\\\ 0&0&0&\cdots&0&\
(\kappa+{\alpha})(\kappa+1+{\alpha})\ &1\end{pmatrix}.$

We can calculate that $\det(\lambda
I-A_{2})=(\lambda-2-{\alpha})(\lambda+1+{\alpha})(\lambda+2+{\alpha})$.
Furthermore, by subtracting $({\kappa}+{\alpha})$ times the $({\kappa}+1)$-th
column from the ${\kappa}$-th column, next pulling out the factor
${\lambda}+{\kappa}+{\alpha}$ from the ${\kappa}$-th column and finally adding
$({\kappa}-1+{\alpha})$ times the ${\kappa}$-th column to the
$({\kappa}-1)$-th column, we can see that $\det(\lambda
I-A_{\kappa})=(\lambda+\kappa+{\alpha})\det(\lambda I-A_{\kappa-1})$.
Therefore, by mathematical induction

$\det(A_{\kappa}-\lambda
I)=\bigl{(}\lambda-(2+{\alpha})\bigr{)}\prod_{l=1}^{\kappa}\bigl{(}\lambda+(l+{\alpha})\bigr{)}.$

We conclude that all eigenvalues are real, and that the only positive
eigenvalue is $\lambda_{1,\kappa}=2+{\alpha}$, so that certainly $\mathop{\rm
Re}\nolimits{\lambda}_{2,{\kappa}}<{\lambda}_{1,{\kappa}}/2$. The transfer
matrix $A_{\kappa}$ is irreducible.

Thus, the conditions of Proposition B.5 are verified. It suffices to identify
the limiting mean.

By (8.6.12) of van der Hofstad (2017), the limiting degree distribution for PA
trees with the PA function $f:k\mapsto k+{\alpha}$ has coordinates

$p_{k}=(2+{\alpha})\frac{\Gamma(k+{\alpha})\Gamma(3+2{\alpha})}{\Gamma(k+3+2{\alpha})\Gamma(1+{\alpha})}.$

This corresponds to the recursion (with $p_{1}=(2+{\alpha})/(3+2{\alpha})$)

$p_{k}=\frac{k-1+{\alpha}}{k+2+2{\alpha}}p_{k-1}.$ (B.9)

We now verify that
$v_{1,\kappa}:=(2+{\alpha})^{-1}\bigl{(}p_{1},p_{2},\dots,p_{\kappa},\sum_{l=\kappa+1}^{\infty}p_{l}(l+{\alpha})\bigr{)}$
is an eigenvector associated with the eigenvalue $2+{\alpha}$ with
$a_{\kappa}^{T}v_{1,\kappa}=1$, and hence obtain that
$n^{-1}X_{n}\rightarrow{\lambda}_{1,{\kappa}}v_{i1,{\kappa}}=\bigl{(}p_{1},p_{2},\dots,p_{\kappa},\sum_{l=\kappa+1}^{\infty}p_{l}(l+{\alpha})\bigr{)}$,
almost surely.

The first coordinate of $A_{\kappa}(2+{\alpha})v_{1,\kappa}$ is
$\sum_{l>1}p_{l}(l+{\alpha})=(2+{\alpha})-p_{1}(1+{\alpha})=(2+{\alpha})p_{1}$,
by Lemma B.3. The second to $\kappa$-th coordinates are equal to
$(2+{\alpha})p_{k}$, for $k=2,\ldots,{\kappa}$, by the relation (B.9). The
$({\kappa}+1)$-st coordinate of $A_{\kappa}(2+{\alpha})v_{1,\kappa}$ is
$\sum_{l>\kappa}p_{l}(l+{\alpha})+(\kappa+{\alpha})(\kappa+1+{\alpha})p_{k}$,
which coincides with $(2+{\alpha})\sum_{l>\kappa}p_{l}(l+{\alpha})$, again by
Lemma B.3.

The covariance function of the limiting Gaussian vector can be obtained from
(B.5) by somewhat tedious calculations, which we omit. We refer to (4.28) of
(Resnick and Samorodnitsky, 2016, page 18) for its exact form.

By the simple linear relation (with drift)
$P_{>k}(n)=1-\sum_{j=1}^{k}P_{j}(n)$, and the continuous mapping theorem, it
follows from the preceding theorem that, for any $\kappa\in\mathbb{N}_{+}$,

$\left(\sqrt{n}(P_{>k}(n)-p_{>k});k=1,\ldots,{\kappa}\right)\rightsquigarrow
N(0,R_{\kappa}),$

The covariance matrix $R_{\kappa}$ in this limit takes a simple form, first
pointed out in Móri (2002), given by

$(R_{\kappa})_{ij}={1}_{\\{i=j\\}}p_{i}(1-p_{i})-{1}_{\\{i\neq
j\\}}p_{i}p_{j}.$ (B.10)

######  Lemma B.3.

The limiting degree distribution $(p_{k})_{k=1}^{\infty}$ in the affine PA
model with the PA function $f(k)=k+{\alpha}$, satisfies, for
$k\in\mathbb{N}_{+}$,

$\sum_{l>k}p_{k}(k+{\alpha})=\frac{(k+{\alpha})(k+1+{\alpha})}{1+{\alpha}}p_{k}.$

###### Proof.

For $k=1$ the left side of the lemma is equal to

$\displaystyle\sum_{l=1}^{\infty}p_{l}(l+{\alpha})-p_{1}(1+{\alpha})$
$\displaystyle=2+{\alpha}-\frac{(2+{\alpha})(1+{\alpha})}{(3+2{\alpha})}=\frac{(1+{\alpha})(2+{\alpha})}{1+{\alpha}}p_{1}.$

This proves the claim for $k=1$. We proceed by mathematical induction. If the
statement is true for any integer up to $k-1$, then the left side of the lemma
is equal to

$\displaystyle\sum_{l>k-1}p_{l}(l+{\alpha})-p_{k}(k+{\alpha})$
$\displaystyle=\frac{(k-1+{\alpha})(k+{\alpha})}{1+{\alpha}}p_{k}\frac{k+2+2{\alpha}}{k-1+{\alpha}}-p_{k}(k+{\alpha}),$

by the induction hypothesis and the relation (B.9) between $p_{k}$ and
$p_{k-1}$. The right side can be reduced to the right side of the lemma.

As a second application of Proposition B.5, we obtain the asymptotic normality
of the empirical degrees in PA models with PA function that is constant
eventually. From our numerical experiments, we infer that the key eigenvalue
condition $\mathop{\rm Re}\nolimits\lambda_{2}(A)<\lambda_{1}(A)/2$ is
generally satisfied when the PA function is sublinear, but we do not know
general conditions for this.

To apply Proposition B.5 to a PA model with eventually constant PA function,
we simply gather all nodes of degrees higher than a cut-off ${\kappa}$ after
which the PA function is constant (not necessarily the smallest such value) in
a single urn, the $({\kappa}+1$)-th one. Since the preferences of the
corresponding nodes are equal, it is not necessary to keep track of the
different degrees of the nodes inside this bin when studying the lower
degrees. The evolution of the empirical degrees
$\bigl{(}P_{1}(n),\ldots,P_{\kappa}(n),P_{>{\kappa}}(n)\bigr{)}$ will be the
same as the evolution of the vectors
$(X_{n,1},\ldots,X_{n,{\kappa}},X_{n,{\kappa}+1})$, if we define the
preferences as $a_{i}=f(i\wedge{\kappa})=f(i)$, and the transition vectors
$\xi_{i}$ for $i\in[{\kappa}]$ by (B.6)–(B.8) with $j$ restricted to
coordinates $j\in[{\kappa}+1]$ and $\xi_{{\kappa}+1}=(1,0,0,\ldots,0)^{T}$.
The last definition corresponds to adding a ball to urn 1 (counting the added
node of degree 1) and moving a ball within the $({\kappa}+1)$-st urn (for
attaching this node to a node of degree $j>{\kappa}$), so not changing
$X_{n,{\kappa}+1}=N_{>{\kappa}}(n)$.

######  Proposition B.4.

In the PA model with PA function $f$ that is constant on $[{\kappa},\infty)$
such that the matrix $A_{\kappa}$ defined in (B.11) satisfies $\mathop{\rm
Re}\nolimits\lambda_{2}(A_{\kappa})<\lambda_{1}(A_{\kappa})/2$, there exist a
probability distribution $(p_{k})$ such that the sequence
$\sqrt{n}\bigl{(}P_{1}(n)-p_{1},\ldots,P_{\kappa}(n)-p_{\kappa},P_{>{\kappa}}(n)-p_{>{\kappa}}\bigr{)}$
tends to a centered normal distribution. If the PA function depends
continuously on a parameter $\theta$, pointwise, then $p_{k}$ and the
covariance matrix depend continuously on ${\theta}$ as well. If the condition
on $A_{\kappa}$ holds for every sufficiently large ${\kappa}$, then the
sequence $\bigl{(}\sqrt{n}\bigl{(}P_{k}(n)-p_{k});k=1,2,\ldots\bigr{)}$
converges in $\mathbb{R}^{\infty}$.

###### Proof.

For the given ${\kappa}$ define an urn process as indicated preceding the
proposition. The matrix $A_{\kappa}$ is given by

$A_{\kappa}=\left(\begin{matrix}0&f(2)&f(3)&f(4)&\cdots&f({\kappa})&f({\kappa}+1)\\\
f(1)&-f(2)\ \,&0&0&\cdots&0&0\\\ 0&f(2)&-f(3)\ \,&0&\cdots&0&0\\\
0&0&f(3)&-f(4)\ \,&\cdots&0&0\\\ 0&0&0&f(4)&\cdots&0&0\\\
\vdots&\vdots&\vdots&\vdots&&\vdots&\vdots\\\ 0&0&0&0&\cdots&-f({\kappa})\
\,&0\\\ 0&0&0&0&\cdots&f({\kappa})&0\\\ \end{matrix}\right).$ (B.11)

The matrix $A_{\kappa}$ is irreducible, and its eigenvalues satisfy
$\mathop{\rm Re}\nolimits\lambda_{2}<\lambda_{1}/2$, by assumption. Thus, the
vector of empirical degrees
$\bigl{(}P_{1}(n),\ldots,P_{\kappa}(n),P_{>{\kappa}}(n)\bigr{)}$, suitably
centered and scaled, is asymptotically normal by Proposition B.5. If this is
true for every sufficiently large ${\kappa}\in\mathbb{N}_{+}$, the infinite
sequence in the final assertion of the proposition converges as well.

It remains to sort out the continuity of the asymptotic mean and covariance if
the PA function depends continuously on a parameter. Denote the parameter by
${\theta}$ and set $q={\kappa}+1$. By its definition the map $\theta\mapsto
A_{\kappa}(\theta)$ inherits the continuity from the PA function. Employing
the $\min$-$\max$ formula (Horn and Johnson (2012, Corollary 8.1.31)),

$\lambda_{1}(\theta)=\max_{x>0}\min_{1\leq i\leq
q}\frac{1}{x_{i}}\sum_{j=1}^{q}\bigl{(}A_{\kappa}(\theta)\bigr{)}_{ij}x_{j},$

where $x>0$ is understood component-wise. By the maximum theorem (e.g. Ok
(2011, page 229)), $\theta\mapsto\lambda_{1}(\theta)$ is continuous. The
corresponding eigenvector can be obtained as
$\operatorname{adj}(C({\theta}))e_{1}$, for $\operatorname{adj}(C)$ the
adjugate matrix of $C=A_{\kappa}-\lambda_{1}I_{q}$. This follows, because
$C\operatorname{adj}(C)=(\det C)I_{q}=0$ and
$\operatorname{adj}(C)e_{1}\not=0$. The latter can be seen from the fact that
the range of $C$, which has dimension $q-1$, is the null space of
$\operatorname{adj}(C)$, as $\operatorname{adj}(C)C=0$, and $e_{1}$ is not in
the range of $C$. Since the adjugate matrix depends continuously on $A$ and
${\lambda}_{1}$, so does the eigenvector
$\operatorname{adj}(C({\theta}))e_{1}$, and this remains valid after scaling.
Inspecting the quantities in (B.2)–(B.4) and (B.5), we see that the asymptotic
covariance matrix is continuous.

##### Acknowledgements

We thank Meiyue Shao for his comments regarding the Perron–Frobenius
eigenvalue of a nonnegative, irreducible matrix.

## References

  * Barabási (2016) Barabási, A.-L. (2016) _Network science_. Cambridge university press.
  * Barabási and Albert (1999) Barabási, A.-L. and Albert, R. (1999) Emergence of scaling in random networks. _science_ , 286, 509–512.
  * Barabási, Albert-László, Albert and Jeong (1999) Barabási, Albert-László, Albert, R. and Jeong, H. (1999) Mean-field theory for scale-free random networks. _Physica A: Statistical Mechanics and its Applications_ , 272, 173–187.
  * Barabási, Albert-László, Albert and Jeong (2000) Barabási, Albert-László, Albert, R. and Jeong, H. (2000) Scale-free characteristics of random networks: the topology of the world-wide web. _Physica A: Statistical Mechanics and its Applications_ , 281, 69–77.
  * Billingsley (2013) Billingsley, P. (2013) _Convergence of probability measures_. John Wiley & Sons.
  * Gao and van der Vaart (2017) Gao, F. and van der Vaart, A. (2017) On the asymptotic normality of estimating the affine preferential attachment network models with random initial degrees. _Stoch. Process. Appl._ , 127, 3754 – 3775.
  * Gao et al. (2017) Gao, F., van der Vaart, A., Castro, R. and van der Hofstad, R. (2017) Consistent estimation in general sublinear preferential attachment trees. _Electron. J. Statist._ , 11, 3979–3999.
  * Hall and Heyde (2014) Hall, P. and Heyde, C. C. (2014) _Martingale limit theory and its application_. Academic press.
  * Horn and Johnson (2012) Horn, R. A. and Johnson, C. R. (2012) _Matrix analysis_. Cambridge university press.
  * Jagers (1975) Jagers, P. (1975) _Branching processes with biological applications_. Wiley.
  * Janson (2004) Janson, S. (2004) Functional limit theorems for multitype branching processes and generalized Pólya urns. _Stochastic Processes and their Applications_ , 110, 177–245.
  * Móri (2002) Móri, T. (2002) On random trees. _Studia Scientiarum Mathematicarum Hungarica_ , 39, 143–155.
  * Nerman (1981) Nerman, O. (1981) On the convergence of supercritical general (C-M-J) branching processes. _Probability Theory and Related Fields_ , 57, 365–395.
  * Ok (2011) Ok, E. A. (2011) _Real analysis with economic applications_. Princeton University Press.
  * Oliveira and Spencer (2005) Oliveira, R. and Spencer, J. (2005) Connectivity transitions in networks with super-linear preferential attachment. _Internet Mathematics_ , 2, 121–163.
  * Resnick and Samorodnitsky (2016) Resnick, S. I. and Samorodnitsky, G. (2016) Asymptotic normality of degree counts in a preferential attachment model. _Advances in Applied Probability_ , 48, 283–299.
  * Rudas, Tóth and Valkó (2007) Rudas, A., Tóth, B. and Valkó, B. (2007) Random trees and general branching processes. _Random Structures & Algorithms_, 31, 186–202.
  * van der Hofstad (2017) van der Hofstad, R. (2017) _Random graphs and complex networks_ , vol. 1. Cambridge University Press.
  * van der Vaart (2000) van der Vaart, A. W. (2000) _Asymptotic statistics (Cambridge series in statistical and probabilistic mathematics)_. Cambridge University Press.
  * van der Vaart and Wellner (1996) van der Vaart, A. W. and Wellner, J. A. (1996) _Weak convergence and empirical processes_. Springer Series in Statistics, Springer-Verlag, New York, with applications to statistics.