Computing the density of the Kesten-Stigum limit in supercritical Galton-Watson processes

TL;DR

This work develops a practical numerical framework to compute the density of the Kesten–Stigum limit $W$ for supercritical Galton–Watson processes with polynomial offspring generating functions. It leverages the Poincaré equation $\varphi(mz)=P(\varphi(z))$ for the Laplace–Stieltjes transform $\varphi(z)=\mathbb{E}[e^{-zW}]$, solving it via forward substitution, a globally convergent fixed-point iteration, a discretized fixed-point scheme, and Newton’s method to obtain moments of $W$. These moments are then used in a moment-matching procedure with a Laguerre–exponential basis, incorporating an atom at $0$ (extinction) and a flexible right tail, with conditioning mitigated through matrix scaling. The approach is validated on several polynomial-offspring examples and applied to real-world ecological data, providing accurate estimates of $W$-density, establishment-time distributions, and generation-size quantiles under non-extinction conditioning. The results offer a scalable, accurate tool for predicting early fluctuations and long-term growth in density-dependent population models and pave the way for extensions to multi-type GW processes.

Abstract

This paper proposes a novel numerical method for computing the density of the limit random variable associated with a supercritical Galton-Watson process. This random variable captures the effect of early demographic fluctuations and determines the random amplitude of long-term exponential population growth. While the existence of a non-trivial limit is ensured by the Kesten-Stigum theorem, computing its density in a stable and efficient manner for arbitrary offspring laws remains a significant challenge. The proposed approach leverages a functional equation that characterizes the Laplace-Stieltjes transform of the limit distribution and combines it with a moment-matching method to obtain accurate approximations within a class of linear combinations of Laguerre polynomials with exponential damping. The effectiveness of the approach is validated on several examples in which the offspring generating function is a polynomial of bounded degree.

Figures (9)

• Figure 1: Relative residual $\mathrm{Res}(\varphi^{(k)})$ along the iterations of the fixed-point method and Newton's method when solving \ref{['eq:poincare']} for $P(z)=\frac{1}{10}(1+z+z^2+\dots+z^9)$.
• Figure 2: Condition numbers as a function of the basis dimension $S$. The blue curve corresponds to $(S+1)\times(S+1)$ Pascal matrices, and the red curve to their row-scaled versions. The yellow curve corresponds to $2(S+1)\times(S+1)$ Pascal matrices, and the purple curve to their scaled counterparts. Row scaling significantly improves conditioning for both square and rectangular matrices. For a fixed basis dimension $S$, the scaled rectangular matrices exhibit the smallest condition numbers and are therefore used in our algorithm
• Figure 3: Approximations of the density of $W$ for GW processes with offspring p.g.f. $P_1(z)$ (left) and $P_2(z)$ (right). In both cases, the extinction probability is $q=0$, and hence there is no point mass at $0$.
• Figure 4: Approximations of the density of $W$ for GW processes with offspring p.g.f. $P_3(z)$ (left) and $P_4(z)$ (right). For $P_3(z)$, a smaller number of generations ($T=8$) was used in the simulation. In both cases, the extinction probability is strictly positive; the red curve represents our approximation of the absolutely continuous part of the density.
• Figure 5: Comparison of the Laguerre-based approximation \ref{['eq:fN']} with the generalized Gamma approximation \ref{['eq:generalized_gamma']} for offspring distributions $P_1(z)$ (left) and $P_3(z)$ (right). Histograms represent empirical densities obtained from simulations.

Theorems & Definitions (13)

• Remark 1
• Theorem 2
• proof
• Lemma 3
• proof
• Remark 4
• Example 5
• Example 6
• Definition 7
• Theorem 8