The field of temporal and contemporaneous aggregations of independent stationary stochastic processes is an important and very active research area in the empirical and theoretical statistics and in other areas as well. The scheme of contemporaneous (also called cross-sectional) aggregation of random-coefficient autoregressive processes of order 1 was firstly proposed by Robinson [16] and Granger [4] in order to obtain the long memory phenomena in aggregated time series. For surveys on papers dealing with the aggregation of different kinds of stochastic processes, see, e.g., Pilipauskaitė and Surgailis [13], Jirak [8, page 512] or the arXiv version of Barczy et al. [2].

In this paper we study the limit behaviour of temporal (time) and contemporaneous (space) aggregations of independent copies of a strictly stationary multitype Galton–Watson branching process with immigration in the so-called iterated and simultaneous cases, respectively. According to our knowledge, the aggregation of general multitype Galton–Watson branching processes with immigration has not been considered in the literature so far. To motivate the fact that the aggregation of branching processes could be an important topic, now we present an interesting and relevant example, where the phenomena of aggregation of this kind of processes may come into play. A usual Integer-valued AutoRegressive (INAR) process of order 1, (Xk)k⩾0 , can be used to model migration, which is quite a big issue nowadays all over the world. More precisely, given a camp, for all k⩾0 , the random variable Xk can be interpreted as the number of migrants to be present in the camp at time k, and every migrant will stay in the camp with probability α∈(0,1) indepedently of each other (i.e., with probability 1−α each migrant leaves the camp) and at any time k⩾1 new migrants may come to the camp. Given several camps in a country, we may suppose that the corresponding INAR processes of order 1 share the same parameter α and they are independent. So, the temporal and contemporaneous aggregations of these INAR processes of order 1 is the total usage of the camps in terms of the number of migrants in the given country in a given time period, and this quantity may be worth studying.

The present paper is organized as follows. In Section 2 we formulate our main results, namely the iterated and simultaneous limit behaviour of time- and space-aggregated independent stationary p-type Galton–Watson branching processes with immigration is described (where p⩾1 ), see Theorems 1 and 2. The limit distributions in these limit theorems coincide, namely, it is a p-dimensional zero mean Brownian motion with a covariance function depending on the expectations and covariances of the offspring and immigration distributions. In the course of the proofs of our results, in Lemma 2.3, we prove that for a subcritical, positively regular multitype Galton–Watson branching process with nontrivial immigration, its unique stationary distribution admits finite αth moments provided that the branching and immigration distributions have finite αth moments, where α∈{1,2,3} . In case of α∈{1,2} , Quine [14] contains this result, however in case of α=3 , we have not found any precise proof in the literature for it, it is something like a folklore, so we decided to write down a detailed proof. As a by-product, we obtain an explicit formula for the third moment in question. Section 3 is devoted to the special case of generalized INAR processes, especially to single-type Galton–Watson branching processes with immigration. All the proofs can be found in Section 4.

Let Z+ , N , R , R+ , and C denote the set of non-negative integers, positive integers, real numbers, non-negative real numbers, and complex numbers, respectively. For all d∈N , the d×d identity matrix is denoted by Id . The standard basis in Rd is denoted by {e1,…,ed} . For v∈Rd , the Euclidean norm is denoted by ‖v‖ , and for A∈Rd×d , the induced matrix norm is denoted by ‖A‖ as well (with a little abuse of notation). All the random variables will be defined on a probability space (Ω,F,P) .

Let (Xk=[Xk,1,…,Xk,p]⊤)k∈Z+ be a p-type Galton–Watson branching process with immigration. For each k,ℓ∈Z+ and i,j∈{1,…,p} , the number of j-type individuals in the kth generation will be denoted by Xk,j , the number of j-type offsprings produced by the ℓth individual belonging to type i of the (k−1)th generation will be denoted by ξk,ℓ(i,j) , and the number of immigrants of type i in the kth generation will be denoted by εk(i) . Then we have (1) Xk= ∑ℓ=1Xk−1,1ξk,ℓ(1,1)⋮ξk,ℓ(1,p)+⋯+ ∑ℓ=1Xk−1,pξk,ℓ(p,1)⋮ξk,ℓ(p,p)+εk(1)⋮εk(p)=: ∑i=1p ∑ℓ=1Xk−1,iξk,ℓ(i)+εk for every k∈N , where we define ∑ℓ=10:=0 . Here {X0,ξk,ℓ(i),εk:k,ℓ∈N,i∈{1,…,p}} are supposed to be independent Z+p -valued random vectors. Note that we do not assume independence among the components of these vectors. Moreover, for all i∈{1,…,p} , {ξ(i),ξk,ℓ(i):k,ℓ∈N} and {ε,εk:k∈N} are supposed to consist of identically distributed random vectors, respectively.

Let us introduce the notations mε:=E(ε)∈R+p , along with Mξ:=E([ξ(1),…,ξ(p)])∈R+p×p and v(i,j):=[Cov(ξ(1,i),ξ(1,j)),…,Cov(ξ(p,i),ξ(p,j)),Cov(ε(i),ε(j))]⊤∈R(p+1)×1 for i,j∈{1,…,p} , provided that the expectations and covariances in question are finite. Let ϱ(Mξ) denote the spectral radius of Mξ , i.e., the maximum of the modulus of the eigenvalues of Mξ . The process (Xk)k∈Z+ is called subcritical, critical or supercritical if ϱ(Mξ) is smaller than 1, equal to 1 or larger than 1, respectively. The matrix Mξ is called primitive if there is a positive integer n∈N such that all the entries of Mξn are positive. The process (Xk)k∈Z+ is called positively regular if Mξ is primitive. In what follows, we suppose that (2) E(ξ(i))∈R+p,i∈{1,…,p},mε∈R+p∖{0},ρ(Mξ)<1,Mξis primitive. For further application, we define the matrix (3) V:=(Vi,j)i,j=1p:=v(i,j)⊤(Ip−Mξ)−1mε1i,j=1p∈Rp×p, provided that the covariances in question are finite.

Under (2), by the Theorem in Quine [14], there is a unique stationary distribution π for (Xk)k∈Z+ . Indeed, under (2), Mξ is irreducible due to the primitivity of Mξ , see Definition 8.5.0 and Theorem 8.5.2 in Horn and Johnson [6]. For the definition of irreducibility, see Horn and Johnson [6, Definitions 6.2.21 and 6.2.22]. Further, Mξ is aperiodic, since this is equivalent to the primitivity of Mξ , see Kesten and Stigum [10, page 314] and Kesten and Stigum [9, Section 3]. For the definition of aperiodicity (also called acyclicity), see, e.g., the Introduction of Danka and Pap [3]. Finally, since mε∈R+p∖{0} , the probability generating function of ε at 0 is less than 1, and E(log( ∑i=1pε(i))1{ε≠0})⩽E( ∑i=1pε(i)1{ε≠0})⩽E( ∑i=1pε(i))= ∑i=1pE(ε(i)), which is finite, so one can apply the Theorem in Quine [14].

For each α∈N , we say that the αth moment of a random vector is finite if all of its mixed moments of order α are finite.

In what follows, we suppose (2) and that the distribution of X0 is the unique stationary distribution π, hence the Markov chain (Xk)k∈Z+ is strictly stationary. Recall that, by (2.1) in Quine and Durham [15], for any measurable function f:Rp→R satisfying E(|f(X0)|)<∞ , we have (4) 1n ∑k=1nf(Xk)⟶a.s.E(f(X0))asn→∞. First we consider a simple aggregation procedure. For each N∈N , consider the stochastic process S(N)=(Sk(N))k∈Z+ given by Sk(N):= ∑j=1N(Xk(j)−E(Xk(j))),k∈Z+, where X(j)=(Xk(j))k∈Z+ , j∈N , is a sequence of independent copies of the strictly stationary p-type Galton–Watson process (Xk)k∈Z+ with immigration. Here we point out that we consider so-called idiosyncratic immigrations, i.e., the immigrations belonging to X(j) , j∈N , are independent.

We will use ⟶Df or Df -lim for weak convergence of finite dimensional distributions, and ⟶D for weak convergence in D(R+,Rp) of stochastic processes with càdlàg sample paths, where D(R+,Rp) denotes the space of Rp -valued càdlàg functions defined on R+ .

We note that using formula (16) presented later on, one could give an explicit formula for Var(X0) (not containing an infinite series).

Note that Propositions 1 and 2 are about the scalings of the space-aggregated process S(N) and the time-aggregated process (∑k=1⌊nt⌋Sk(1))t∈R+ , respectively.

For each N,n∈N , consider the stochastic process S(N,n)=(St(N,n))t∈R+ given by St(N,n):= ∑j=1N ∑k=1⌊nt⌋(Xk(j)−E(Xk(j))),t∈R+.

A key ingredient of the proofs is the fact that (Xk−E(Xk))k∈Z+ can be rewritten as a stable first order vector autoregressive process with coefficient matrix Mξ and with heteroscedastic innovations, see (24).

We devote this section to the analysis of aggregation of Generalized Integer-Valued Autoregressive processes of order p∈N (GINAR(p) processes), which are special cases of p-type Galton–Watson branching processes with immigration introduced in (1). For historical fidelity, we note that it was Latour [11] who introduced GINAR(p) processes as generalizations of INAR(p) processes. This class of processes became popular in modelling integer-valued time series data such as the daily number of claims at an insurance company. In fact, a GINAR(1) process is a (general) single type Galton–Watson branching process with immigration.

Let (Zk)k⩾−p+1 be a GINAR(p) process. Namely, for each k,ℓ∈Z+ and i∈{1,…,p} , the number of individuals in the kth generation will be denoted by Zk , the number of offsprings produced by the ℓth individual belonging to the (k−i)th generation will be denoted by ξk,ℓ(i,1) , and the number of immigrants in the kth generation will be denoted by εk(1) . Here the 1-s in the supercripts of ξk,ℓ(i,1) and εk(1) are displayed in order to have a better comparison with (1). Then we have Zk= ∑ℓ=1Zk−1ξk,ℓ(1,1)+⋯+ ∑ℓ=1Zk−pξk,ℓ(p,1)+εk(1),k∈N. Here {Z0,Z−1,…,Z−p+1,ξk,ℓ(i,1),εk(1):k,ℓ∈N,i∈{1,…,p}} are supposed to be independent nonnegative integer-valued random variables. Moreover, for all i∈{1,…,p} , {ξ(i,1),ξk,ℓ(i,1):k,ℓ∈N} and {ε(1),εk(1):k∈N} are supposed to consist of identically distributed random variables, respectively.

A GINAR(p) process can be embedded in a p-type Galton–Watson branching process with immigration (Xk=[Zk,…,Zk−p+1]⊤)k∈Z+ with the corresponding p-dimensional random vectors ξk,ℓ(1)=ξk,ℓ(1,1)10⋮0,…,ξk,ℓ(p−1)=ξk,ℓ(p−1,1)0⋮01,ξk,ℓ(p)=ξk,ℓ(p,1)00⋮0,εk=εk(1)00⋮0 for any k,ℓ∈N .

In what follows, we reformulate the classification of GINAR(p) processes in terms of the expectations of the offspring distributions.

Next, we specialize the matrix V , defined in (3), in case of a subcritical GINAR(p) process.

Finally, we specialize the limit distribution in Theorems 1 and 2 in case of a subcritical GINAR(p) process. Remark 5.

Let us note that in case of p=1 and E(ξ(1,1))<1 (yielding that the corresponding GINAR(1) process is subcritical), the limit process in Theorems 1 and 2 can be written as 11−E(ξ(1,1))E(ε(1))Var(ξ(1,1))+(1−E(ξ(1,1)))Var(ε(1))1−E(ξ(1,1))W, where W=(Wt)t∈R+ is a standard 1-dimensional Brownian motion. Indeed, this holds, since in this special case Mξ=E(ξ(1,1)) yielding that (Ip−Mξ)−1=(1−E(ξ(1,1)))−1 , and, by (10), V=Cov(ξ(1,1),ξ(1,1))Cov(ε(1),ε(1))⊤E(ε(1))1−E(ξ(1,1))1=Var(ξ(1,1))E(ε(1))1−E(ξ(1,1))+Var(ε(1)).□