Population Dynamics
You’ve very likely seen a population pyramid before. I thought it could be interesting to write down a mathematical description of the time evolution of such diagrams, and simulate it.
This post starts by showing the PDE governing the evolution of the population, and how to discretize it & various probability distributions required. Then follows an interactive simulation of a population over time, using the preceding scheme. Lastly, I attempt to mathematically analyze the equation, to see under what conditions various behaviours arise.
Contents
Mathematical Framework
Before diving in, I will make a couple of simplifying assumptions:
- The male & female populations are always balanced; no need to track them separately.
- There is no immigration; the population only grows by birth.
- Fractional people exist; the populations are given by real numbers, not integers.
These assumptions could easily be relaxed, but I don’t think that would make the problem more interesting mathematically.
Population PDE
We would like to model the ‘demographic curve’ - \(n(a, t)\) - which tells us how many people of any given age - \(a\) - are alive at time \(t\).
The time evolution of \(n\) has two contributing terms. Firstly an aging term; everybody gets older at the same rate, namely one year per year. Meaning we will have a right-moving advection term with velocity equal to 1. Secondly a hazard term, representing mortality.
$$\boxed{ \pfrac{n}{t} = -\pfrac{n}{a} - h(a)n \label{eqn_dndt}\tag{1} }$$The hazard rate - \(h(a)\) - is the probability density per unit time that a person of a given age dies (or emigrates, perhaps).
Lastly we need a boundary condition for \(a=0,\) because we cannot have a negative age to advect from. The boundary condition is imposed by the birth-rate
$$n(0, t) = B(t) \label{eqn_n_bc}\tag{2} $$We must now define \(B(t).\)
Birth & Death Rates
The endogenous birth-rate (density) - \(b(a)\) - is the probability density per unit time that a person of a given age reproduces. Therefore the total number of endogenous births per unit time is given by
$$B_n(t) = \int_0^\infty b(a) n(a, t) \,da \label{eqn_B}\tag{3} $$The exogenous birth-rate - \(B_x(t)\) - is a birth-rate which is independent of the population. The total number of births per unit time is their sum \(B = B_n + B_x.\)
Similarly, the total number of deaths per unit time is
$$D(t) = \int_0^\infty h(a) n(a, t) \,da \label{eqn_D}\tag{4} $$Total Population
The total population is given by
$$P(t) = \int_0^\infty n(a, t) \,da$$What is the rate of change of the total population? Well …
$$ \begin{flalign} && \frac{dP}{dt} &= \frac{d}{dt}\int_0^\infty n \,da & \text{Definition of }p \\ && &= \int_0^\infty \pfrac{n}{t}\,da & \text{Linearity} \\ && &= -\int_0^\infty \Big(\pfrac{n}{a} + n h\Big) \,da & \text{Equation }\eqref{eqn_dndt} \\ && &= n(0, t) - n(\infty, t) - \int_0^\infty nh\,da & \text{Fundamental Thm.}\vphantom{\pfrac{}{}} \\ && &= B - D & \text{Equations \eqref{eqn_D} and \eqref{eqn_n_bc}} \vphantom{\pfrac{}{}} \end{flalign} $$To go from the penultimate to the final line we have also assumed that there are no infinitely old people. Feels like a reasonable assumption. Thus the rate of change of the total population is the total births minus the total deaths. Unsurprising but a good check that things are thus far consistent.
Survival, Mortality & Hazard Curves
Equation \eqref{eqn_dndt} is written in terms of the hazard rate, which is a probability density. However, for our discretized simulation, we need actual probabilities and not densities. In this section we will derive the conversion from the hazard rate to actual probabilities, and relate them both to what fraction of people survive to each age.
The contents of this section are related to survival analysis.
Let’s start by defining the mortality function - \(M(a).\) This is the probability that an individual has died by age \(a.\) It is the inverse of the survival function \(S(a) = 1 - M(a),\) which is the probability that an individual has not died by age \(a.\)
Say a person lives to age \(A,\) then \(M(a) = \pp (A \leq a).\) The unconditional (a priori) probability density of dying at age \(a\) is
$$m(a) = \dfrac{M}{a}$$This has the property that it is small for large ages, which is perhaps confusing at first sight but makes perfect sense. What we care about is not the a priori probability density of dying at some age, but rather the probability density conditioned on having survived up to that point. This quantity is our hazard rate.
Intuitively, we can see that if we survived to age \(a_1,\) then the probability of dying between ages \(a_1\) and \(a_2 \geq a_1\) is a simple rescaling of our unconditional mortality curve:
$$H(a_2|a_1) = \frac{M(a_2) - M(a_1)}{1 \,\, - \,\, M(a_1)} = \frac{S(a_1) - S(a_2)}{S(a_1)} \label{eqn_H}\tag{5} $$Thus \(H(a_2|a_1)\) is the hazard probability.
Example
Consider a simple mortality function given by
$$M(a) = \mathrm{erf}\big((a - 80) / 5\big)$$This would tell us, for example, that the a priori probability density of dying at age 100 is quite low; \(M(101) - M(100) \approx 1.8\times 10^{-5}.\) This does not mean 100 year olds are all hale and sprightly, it means we were simply very unlikely to survive to 100 in the first place. The hazard of being aged 100 is actually quite high!
$$H(101|100) = \frac{M(101) - M(100)}{1 - M(100)} \approx \frac{1.8\times 10^{-5}}{3.2 \times 10^{-5}} = 0.57$$The hazard rate is given by the derivative of the hazard probability wrt \(a_2\) evaluated at \(a_1:\)
$$h(a_1) = \pfrac{H(a_2|a_1)}{a_2} \bigg|_{a_2 = a_1} = -\frac{1}{S(a_1)}\dfrac{S(a_2)}{a_2}\bigg|_{a_2 = a_1} \label{eqn_h_rate}\tag{6} $$i.e.
$$\dfrac{S}{a} = -h(a)S$$This is a simple ODE which has solution
$$S(a) = C\exp\left({-\!\int_0^a h(a^\prime) \,da^\prime}\right)$$We can pin down the value of the integration constant \(C\) by noting that \(S(0)\equiv 1\) thus \(C = 1,\) so
$$S(a) = \exp\left({-\!\int_0^a h(a^\prime) \,da^\prime}\right) \label{eqn_S_cont}\tag{7} $$Thus, we can either assume that the survival function is given and use equation \eqref{eqn_h_rate} to calculate the hazard rate. Or, vice-versa, we can assume we are given the hazard rate and use equation \eqref{eqn_S_cont} to calculate the survival function. It doesn’t really matter.
Life Expectancy
The life-expectancy is the mean age at which death occurs i.e.
$$\mu_A = \mathbb E[A] = \int_0^\infty a m(a) \,da$$Solving the Equations Numerically
Discretization & Update Scheme
For our simulation we have both a discrete grid for our ages, and discrete time steps. First, lets choose a grid of ages to represent, with step size \(\Delta a\). Thus we say we know the value of \(n\) at ages \(a_i = i\Delta a\) for \(i \geq 0.\) We require the population to age \(\Delta a\) each time step, so that it lands exactly on our grid, and we suffer no numerical diffusion. Therefore the grid spacing \(\Delta a\) also sets the time step and we have \(t_i = i\Delta a.\)
In the time discretized setting, we cannot work with the hazard rate, but rather we work with the previously defined hazard probability. Specifically we need the probability of dying in any given time step given we reached a certain age. So, rather than considering \(h(a)\) we use \(H(a+\Delta a|a).\) In this way we turn the 2d function \(H\) into the 1d function of \(a\) needed for our simulation. Further, define \(H_i = H(a_i+\Delta a | a_i)\) to be the value of this function at grid point \(i.\)
In our simulation it is entirely possible to have a probability of 1 of dying in the next \(\Delta a\) time. This is also possible in the continuous setting but it implies an infinite hazard rate at the point the survival curve becomes zero. You can see this either by noting \(0 = e^{-\infty}\) from equation \eqref{eqn_S_cont} or the fact we get zero denominator in equation \eqref{eqn_H}
The endogenous birth-rate \(b\) can also be discretized simply by integrating over timesteps;
$$b_i = \int_{a_i}^{a_i+\Delta a} b(a) \,da$$The total birth-rate then becomes a sum
$$B(t_j) = \sum_{i=0}^\infty b_i n(a_i, t_j)$$Our final discretized update scheme is very simple
$$\begin{align} n(a_i, t_j) &= n(a_{i-1}, t_{j-1}) H_i & \qquad i > 0\\ n(a_0, t_j) &= B(t_{j}) H_0 & \end{align}$$but it exactly captures the aging (advection) term and the hazard term of the continuous setting.
Discrete Survival Curve
As mentioned above, we have a choice; we can consider the hazard probabilities to be given, and use these to calculate the survival function, or we can calculate the hazard probabilities from the survival function.
If we consider the hazard probabilities \(H_i\) to be given, we can use them to calculate \(S_i\) thus:
$$S_i = \prod_{j=0}^{i-1} (1 - H_i)$$and we already know that if we consider \(S_i\) to be given then we calculate \(H_i\) as
$$H_i = \frac{S_{i+1} - S_i}{S_i}$$These equations actually hold for any partition of the interval \([0, a_i]\) not just our regular grid.
Interactive Simulation
We’re ready to assemble this together into a simulation of our population evolution. Use the slider inputs below the simulation to play with the birth- & hazard rates, and see what happens. The code for this simulation is available on GitHub.
Analysis of the Equations & Population Evolution
A very natural question to ask is how does the population behave? Do there exist solutions where \(n(a, t)\) is not changing with time? What about the situation where \(p(t)\) is not changing with time? If such solutions exist, are they stable?
Reproduction Rate
One obvious criterion for the solution to be steady-state is that every person should produce (on average) exactly one offspring over their lifetime. Thereby replacing themself, and their child then replacing themself, and so on ad-infinitum.
We know that the probability density of reproducing is \(b(a)\), so we may initially think the overall average number of children per individual is the integral of this quantity. It almost is. But a person may have already succumbed by age \(a,\) so the birth-rate density at this age “counts less”. How much less? It gets weighted by the probability of surviving (at least) to age \(a\) - \(S(a)\).
The integral of this weighted birth-rate gives us the (endogenous) reproduction rate.
$$R = \int_0^\infty b(a) S(a) \,da \label{eqn_R}\tag{8} $$Note, for this value to be finite we require that \(b(a)S(a) \rightarrow 0\) fast enough that the integral converges.
Parameter Cases
Remembering that we also have to consider the effect of the exogenous birth-rate we get 5 cases to consider.
| Case | Parameters | Behaviour |
|---|---|---|
| 1 | \(\quad R > 1 \quad\) | \(P\) will grow without bound. |
| 2 | \(\quad R = 1, B_x > 0 \quad\) | \(P\) will grow without bound. |
| 3 | \(\quad R = 1, B_x = 0 \quad\) | \(P\) is stable\(^{[1]}\) (but not necessarily steady-state). |
| 4 | \(\quad R < 1, B_x > 0 \quad\) | \(P\) will approach a steady state. |
| 5 | \(\quad R < 1, B_x = 0 \quad\) | \(P\) will shrink to zero over time.\(^{[2]}\) |
\(^{[1]}\) It neither blows-up, nor shrinks to zero but can be “oscillating” around a long-run value.
\(^{[2]}\) Thus it approaches the trivial steady state solution of no population!
For any value of \(R\), if \(B_x=0\) then \(n(a, t) = 0\) is a valid solution. This trivial case is uninteresting, I merely note that it exists.
By playing with the interactive simulation above you can see all 5 of these cases display the listed behaviour. The horizontal line marked ‘1’ on the right-hand y-axis corresponds to \(R=1\) for the birth-rate. Next we can try to prove if we must see these behaviours, or under which conditions they arise.
Lagrangian Formulation
We can view equation \eqref{eqn_dndt} through the lense of the total derivative. For a cohort born in year \(\tau\), their age changes as \(a(t) = t - \tau.\) Thus
$$\frac{d n}{d a} = \pfrac{n}{a} + \pfrac{t}{a}\pfrac{n}{t} = \pfrac{n}{a} + \pfrac{n}{t} = - h(a)\cdot n \label{eqn_dnda_tot}\tag{9} $$This tells us that, for a specified cohort, the population evolution is an governed by nothing more than an exponential decay with the hazard rate. This looks familiar; it is the same differential equation we saw for the survival curve \(S.\) Thus the solution is
$$n = C \exp\left(-\int_{0}^a h(a^\prime)\,da^\prime\right) = C S(a)$$We know the initial population of this cohort was given by the birth-rate at time \(\tau\), which sets the integration constant \(C = B(\tau).\) Thus
$$n(a, t) = B(t - a) S(a) \label{eqn_n_soln_1}\tag{10} $$This is a useful equation because it lets us rewrite various quantities involving \(n\) in terms of \(B.\) For example
$$P(t) = \int_0^\infty n(a, t) \,da = \int_0^\infty B(t-a)S(a) \,da$$This equation makes intuitive sense, the total population now is given by contributions from the birth-rate at each previous time, weighted by how many people survived form those previous times until now.
Further, combining equation \eqref{eqn_n_soln_1} with equation \eqref{eqn_B} gives us an equation solely in terms of the birth-rate:
$$B_n(t) = \int_0^\infty b(a)S(a) B(t-a)\,da \label{eqn_B_soln_1}\tag{11} $$Aside
I had somehow expected we would be able to write things as a delay differential equation for \(B,\) however that isn’t quite what we get. You can take the time derivative of equation \eqref{eqn_B_soln_1} and you get
$$\dfrac{B_n}{t} = \int_0^\infty b(a)S(a) \pfrac{B}{t}(t-a)\,da$$So the rate of change of \(B_n\) (and thus \(B\)) depends not on its values in the past, but rather on how it was changing in the past.
In any case, equation \eqref{eqn_B_soln_1} looks almost like the definition of \(R\) (see equation \eqref{eqn_R}), which makes it useful for characterizing the behaviour.
Existence of Steady-State Solutions
Equation \eqref{eqn_B_soln_1} actually lets us directly show for which cases a steady-state solution does exist. At steady state we must have constant endogenous and exogenous birth-rates. Because the birth-rates are constant then we can pull them out of the integral in \eqref{eqn_B_soln_1} giving
$$B_n = (B_n + B_x)R$$If \(B_x = 0,\) this means we have \(B_n = R B_n,\) which is true only when \(R = 1.\) Thus for \(R = 1\) and \(B_x = 0\) we have a steady-state solution, this is case 3 above.
Now consider \(B_x > 0.\) Rearranging slightly we get
$$B_n = B_x \frac{R}{1-R}$$This clearly only has a physically meaningful solution if \(R < 1,\) which corresponds to case 4 above. Each person does not quite replace themselves, but this is compensated for by a non-zero exogenous birth-rate. As \(R\) approaches one from below, \(B_n\) blows up towards infinity. This means the steady state total population gets very large as our reproduction rate approaches unity.
For \(R = 1\) the right-hand side is undefined and there is no steady-state solution. The reason is clear; at \(R = 1\) each person exactly replaces themselves via the endogenous birth-rate, but we are also constantly increasing the population via the exogenous birth-rate.
For \(R > 1\) the endogenous birth-rate becomes negative which is physically meaningless. There is no constant \(B_n > 0\) which solves equation equation \eqref{eqn_B_soln_1} in this case.
Note, this says only that steady-state solutions exist in the two cases specified; it does not say that we will always approach these solutions from any set of initial conditions.
For example, in the case of \(B_x = 0, R=1,\) if \(b\) is given by a delta function, then we will get delta-pulse baby-booms, which persist forever. However we will never approach a true steady-state solution. I believe, on the other hand, that if there is some width to \(b\) then we will approach the steady-state solution but I have not proved that.
Expanding & Contracting Total Population
For the cases without a steady-state solution (cases 1, 2, and 5) it seems difficult to make any concrete statements with full generality.
A delta pulse of birth-rate (i.e. a single cohort) will be spread by convolution with the endogenous birth-rate \(b,\) which itself gets convolved with \(b,\) spreading it further etc. It seems intuitive that, for \(R > 1,\) eventually there will have been more people born than contained in the initial cohort. However, the convolution means they can be spread out so that at any one moment the birth-rate is lower than the initial rate.
In fact, its even possible for the population to be non-monotonic but still grow over time; consider a situation where 20% of people die in infancy, but then adults have on average 2 children (per parent). In this situation each person has \(0.8 \times 2 = 1.6\) offspring, and the population clearly grows with each generation, despite initially falling due to the infant mortality.
Probably because of such “corner cases”, I haven’t found a convincingly rigorous argument for the “facts” that, in the long run, the population must grow if \(R > 1\) and shrink if \(R < 1\). I suspect that any proof will rely on the fact that there is not infinite “space” for the generational convolutions to expand into. The survival curve allows for a certain range of ages, so eventually the trailing edge of one generation starts overlapping with leading edge of the next. This means we probably must specify some properties that \(S(a)\) and \(S(a)b(a)\) obey, in order to obtain a proof.
Cohort Resultant Population
A last note.
Effectively what we are asking involves counting the current (live) population arising from a single cohort at each moment in time. That is to say, how many of the cohort itself are alive, how many of the cohorts direct offspring, their offspring, and so forth? Let’s write that down, at least.
Without loss of generality, consider a cohort size of 1; due to linearity of integrals any other size is just given by a factor on the whole expression. For a cohort of age \(a\) the total current population arising due to that cohort is given by
$$ p(a) = S(a) + \int_0^{a} b(a^\prime) S(a^\prime) p(a - a^\prime)\,da^\prime\qquad t \geq \tau \label{eqn_p_cohort}\tag{12} $$This expression can be easily understood. The first term is how many of the original cohort are surviving at age \(a.\) The second is the current resultant population arising from each direct child cohort, whose age is the parent cohort’s current age, less the parents’ age when the children were born. Again, \(p\) represents a “unit” cohort and it is weighted by the birth-rate when they were born. The recursive nature of the expression captures subsequent generations arising from the original cohort.
What we are seeking to prove is that \(p(a)\) diverges for \(R > 1\) and converges to zero for \(R < 1.\)
In principle, equation \eqref{eqn_p_cohort} really tells us a lot about the dynamics of the total population; it captures everything about the behaviour that arises due to the endogenous birth-rate. If we could solve it, the whole dynamics of the total population - \(P\) - could be obtained by convolving \(p\) with the exogenous birth-rate \(B_x\). (This is not quite as much information as is contained in the demographic curve - \(n\) - of course, but still good.)