Stable Populations
In this unit we do some stable population calculations, including determination of the intrinsic growth rate, as illustrated in Box 7.1, and computing he stable equivalent age distribution, as illustrated in Box 7.2 in the textbook.
I used Mata throughout, which I think leads to clearer code,
but I need it only for the eigenanalysis in the third section;
otherwise all calculations can be done in plain Stata,
as shown in the alternative version.
Box 7.1. The Population of Egypt
We showed in the previous unit how to calculate r from the first eigenvalue of the Leslie Matrix. We now use the Egyptian example in Box 7.1 to illustrate Coale's method. We start by entering person-years lived for ages 15-19 to 45-49, the maternity function at those ages, and the midpoints of the age groups.
. mata: ------------------------------------------------- mata (type end to exit) ------------------------------------------- : L = (4.66740, 4.63097, 4.58518, 4.53206, 4.46912, 4.39135, 4.28969) : m = (0.00567, 0.06627, 0.11204, 0.07889, 0.05075, 0.01590, 0.00610) : a = (15,20,25,30,35,40,45) :+ 2.5 : end ---------------------------------------------------------------------------------------------------------------------
The Net Reproduction Ratio
The Net Reproduction Ratio NRR is easily computed as the sum of the products of the survival ratios and the maternity function
. mata: ------------------------------------------------- mata (type end to exit) ------------------------------------------- : nrr = sum( L :* m ) : nrr 1.527413734 : end ---------------------------------------------------------------------------------------------------------------------
The NRR is 1.527 daughters per woman, in agreement with the textbook.
Coale's Method for Estimating r
Next we solve Lotka's equation. Coale's method follows from a Taylor expansion of Lotka's integral (equation 7.10b in the text) which I'll call f(r), and its first derivative f'(r)=-f(r)A(r) where A is the mean age of childbearing written as a function of r. If r0 is a trial value and r the true value or solution, write f(r0) = f(r) + (r0-r)f'(r) and note that f(r)=1 and f'(r)=-A where A is the true mean age of childbearing in the stable population. The next trial value would then be r = r0 + (f(r0)-1)/A. As the textbook notes, we don't know A until we know r, so we use an approximation, in this example 27. For a starting value we use r = log(NRR)/A.
We write a one-line function to compute a discrete approximation to Lotka's integral given a trial value of r. We then use log(NRR)/27 as a trial value and compare the result to unity to determine a new trial value. The code below does three iterations
. mata:
------------------------------------------------- mata (type end to exit) -------------------------------------------
: scalar f(scalar r, vector a, vector L, vector m) {
> return( sum(exp(-r * a) :* L :* m ) )
> }
: // starting value
: r = log(nrr)/27
: fr = f(r,a,L,m)
: r,fr
1 2
+-----------------------------+
1 | .0156879976 .9583727822 |
+-----------------------------+
: // three iterations
: for(niter=1; niter<=3; niter++) {
> r = r + (f(r,a,L,m)-1)/27
> fr = f(r,a,L,m)
> r,fr
> }
1 2
+-----------------------------+
1 | .0141462488 1.002887 |
+-----------------------------+
1 2
+-----------------------------+
1 | .0142531747 .9997313309 |
+-----------------------------+
1 2
+-----------------------------+
1 | .014243224 1.000024566 |
+-----------------------------+
: end
---------------------------------------------------------------------------------------------------------------------
The three iterations are the same as in the textbook. The intrinsic growth rate of Egypt in 1997 was .01425.
Using Exact First Derivatives
It occurred to me studing Coale's method that one could have expanded Lotka's integral around the trial value instead of the true value. Write f(r)=f(r0) + (r-r0)f'(r0) = f(r0) - (r-r0)f(r0)A(r0). This requires evaluating the mean age of childbearing at the trial value, which we can always do, and leads to Newton's method, with improved trial value r = r0 + (f(r0)-1)/(f(r0)A(r0)). The procedure converges in one iteration in this example:
. mata:
------------------------------------------------- mata (type end to exit) -------------------------------------------
: scalar g(scalar r, vector a, vector L, vector m) {
> return( sum(a :* exp(-r * a) :* L :* m ) )
> }
: // starting value
: log(nrr)/27
.0156879976
: fr = f(r,a,L,m)
: gr = g(r,a,L,m)
: r,fr
1 2
+-----------------------------+
1 | .014243224 1.000024566 |
+-----------------------------+
: // first iteration
: A = gr/fr
: r = r + (fr-1)/(fr*A)
: fr = f(r,a,L,m)
: r,fr,g(r,a,L,m)/fr
1 2 3
+-------------------------------------------+
1 | .0142440575 1 29.47248278 |
+-------------------------------------------+
: end
---------------------------------------------------------------------------------------------------------------------
As a bonus we get the mean age of childbearing in the stable population, which is 29.47.
Lotka used a quadratic expansion to obtain a non-iterative approximate solution.
Box 7.2 The Population of the U.S.
Let us now turn to the example in Box 7.2, which involves the female
population of the U.S. in 1991. The data include the actual age
distribution as well as the person-years lived and the maternity
function. These are available in a file called
prestonBox72.dat in the datasets section.
. infile age ca L m using ///
> http://data.princeton.edu/eco572/datasets/prestonBox72.dat, clear
(18 observations read)
. list
+-------------------------------+
| age ca L m |
|-------------------------------|
1. | 0 .0726 4.95804 0 |
2. | 5 .0689 4.95002 0 |
3. | 10 .0667 4.94603 .0007 |
4. | 15 .0648 4.93804 .0303 |
5. | 20 .0729 4.92552 .0566 |
|-------------------------------|
6. | 25 .0799 4.91138 .0578 |
7. | 30 .0861 4.89356 .0388 |
8. | 35 .0801 4.86941 .0157 |
9. | 40 .0735 4.83577 .0027 |
10. | 45 .0556 4.78475 .0001 |
|-------------------------------|
11. | 50 .0464 4.70374 0 |
12. | 55 .0421 4.57712 0 |
13. | 60 .0436 4.38502 0 |
14. | 65 .0429 4.10756 0 |
15. | 70 .0365 3.7199 0 |
|-------------------------------|
16. | 75 .0294 3.19192 0 |
17. | 80 .0203 2.49203 0 |
18. | 85 .0176 2.73044 0 |
+-------------------------------+
Eigen-Analysis of the Leslie Matrix
The first step will be to construct a Leslie matrix and carry out the eigenanalysis. This is easy to do using the function we wrote in the last unit. (If you are following along in a different session you can read the function from a do file, as I do below.)
. quietly do http://data.princeton.edu/eco572/leslie.do . mata: ------------------------------------------------- mata (type end to exit) ------------------------------------------- : L = st_data(.,"L") : m = st_data(.,"m") : M = Leslie(L,m) : values = J(0,0,.) : vectors = J(0,0,.) : eigensystem(M,vectors,values) : values[1] .99916795 : log(values[1])/5 -.000166479 : sa = Re(vectors[,1]/sum(vectors[,1])) : end ---------------------------------------------------------------------------------------------------------------------
The intrinsic rate of growth of the U.S. in 1991 was r = -0.0001665. If the observed fertility and mortality patterns were to prevail the population of the U.S. would (eventually) decline about 0.017% per year.
Estimating r using Coale's Method
Let us now calculate the rate using Coale's method to solve Lotka's equation iteratively. We start by computing the NRR, which as you would expect is less than one. We also need the midpoints of the age groups. For the last group we use 6.79, life expectancy at age 85, as in the textbook
. mata:
------------------------------------------------- mata (type end to exit) -------------------------------------------
: nrr = sum( m :* L )
: nrr
.995601948
: r = log(nrr)/27
: a = range(0,85,5) :+ 2.5
: a[18] = 85 + 6.79
: fr = f(r,a,L,m)
: r,fr
1 2
+-------------------------------+
1 | -.0001632501 .9999144141 |
+-------------------------------+
: // three iterations
: for(niter=1; niter<=3; niter++) {
> r = r + (f(r,a,L,m)-1)/27
> fr = f(r,a,L,m)
> r,fr
> }
1 2
+-------------------------------+
1 | -.0001664199 .9999983443 |
+-------------------------------+
1 2
+-------------------------------+
1 | -.0001664812 .999999968 |
+-------------------------------+
1 2
+-------------------------------+
1 | -.0001664824 .9999999994 |
+-------------------------------+
: end
---------------------------------------------------------------------------------------------------------------------
As you can see the method settles on r = -0.0001665, in agreement with the earlier result. Let us try using the exact first derivative
. mata:
------------------------------------------------- mata (type end to exit) -------------------------------------------
: r = log(nrr)/27
: fr = f(r,a,L,m)
: gr = g(r,a,L,m)
: r,fr,gr/fr
1 2 3
+----------------------------------------------+
1 | -.0001632501 .9999144141 26.47876384 |
+----------------------------------------------+
: // one iteration
: A = gr/fr
: r = r + (fr-1)/(fr*A)
: fr = f(r,a,L,m)
: gr = g(r,a,L,m)
: r,fr,gr/fr
1 2 3
+----------------------------------------------+
1 | -.0001664826 1.000000004 26.4788846 |
+----------------------------------------------+
: end
---------------------------------------------------------------------------------------------------------------------
The method converges in a single iteration (another crank would change the result by one in the tenth decimal place), and gives a mean age of childbearing in the stable population of 26.479. Note that there is a typo in the textbook, which reports r as -.00018 on page 150.
Obtaining the Equivalent Age Distribution
The final step will be to calculate the stable equivalent age distribution. This is a simple function of r and the person-years lived. (It also involves the birth rate, but that is also a function of r and person-years lived.) Below we compute b and then the age distribution, and store the two calculations in Stata for comparison and plotting
. mata:
------------------------------------------------- mata (type end to exit) -------------------------------------------
: b = 1/sum( exp(-r*a) :* L )
: b
.012584262
: saa = b * exp(-r*a) :* L
: st_store(.,st_addvar("float","sa"), sa)
: st_store(.,st_addvar("float","saa"), saa)
: end
---------------------------------------------------------------------------------------------------------------------
. list age ca sa saa
+-----------------------------------+
| age ca sa saa |
|-----------------------------------|
1. | 0 .0726 .0624188 .0624192 |
2. | 5 .0689 .0623698 .0623702 |
3. | 10 .0667 .0623714 .0623718 |
4. | 15 .0648 .0623225 .0623229 |
5. | 20 .0729 .0622162 .0622167 |
|-----------------------------------|
6. | 25 .0799 .0620893 .0620897 |
7. | 30 .0861 .0619155 .0619159 |
8. | 35 .0801 .0616613 .0616617 |
9. | 40 .0735 .0612863 .0612867 |
10. | 45 .0556 .0606902 .0606906 |
|-----------------------------------|
11. | 50 .0464 .0597123 .0597127 |
12. | 55 .0421 .0581533 .0581537 |
13. | 60 .0436 .055759 .0557594 |
14. | 65 .0429 .0522744 .0522748 |
15. | 70 .0365 .0473803 .0473806 |
|-----------------------------------|
16. | 75 .0294 .0406893 .0406896 |
17. | 80 .0203 .0317938 .0317941 |
18. | 85 .0176 .0348964 .0348897 |
+-----------------------------------+
As you can see the do calculations of the stable equivalent age distribution give essentially the same result. Watcher's textbook explains why this is so.
Let us now plot the current and stable age distributions using the midpoints of the age groups
. mata st_store(.,st_addvar("float","am"), a)
. line ca sa am, lpat(dash solid) ///
> xtitle("age") ytitle(Prop in 5-year age group) ///
> title(Actual and Stable Equivalent Age Distributions) ///
> subtitle(U.S. females 1991) ///
> legend(order(1 "Actual" 2 "Stable") col(1) ring(0) pos(1))
. graph export stablepop.png, replace
(file stablepop.png written in PNG format)
The graph reproduces Figure 7.4 in the textbook, showing how the current age distribution has proportionately more people at younger ages (up to 40-44) and fewer at older ages (45-49 and up).
A Preview of Momentum
The 'youthful' U.S. age structure means that, at current rates, the U.S. population will continue to grow for a while before it starts to decline, so it has positive population momentum. How long would it grow and what size would it reach? The following code snipet answers both questions.
. mata:
------------------------------------------------- mata (type end to exit) -------------------------------------------
: delta = 1
: c0 = st_data(.,"ca") // current age distribution
: n=0
: first = 1
: while (delta > 1e-8) {
> c1 = M * c0
> r1 = log(sum(c1)/sum(c0))/5
> delta = abs(r1-r)
> n++
> // print when pop first declines
> if(r1 < 0 & first) {
> n, sum(c1), r1, delta
> first = 0
> }
> c0 = c1
> }
1 2 3 4
+-------------------------------------------------------------+
1 | 9 1.138400994 -.0002235473 .0000570647 |
+-------------------------------------------------------------+
: n, sum(c1), r1
1 2 3
+----------------------------------------------+
1 | 53 1.084463515 -.0001664865 |
+----------------------------------------------+
: end
---------------------------------------------------------------------------------------------------------------------
We see that the population grows for 9 projection periods (45 years), reaching a size 13.84% larger than today, and then heads for extinction, with the growth rate reaching -0.0001665 (to a close approx) after 56 projection periods (280 years), at which time it is still 8.45% larger than at the start. (So extinction is not imminent!)
Question: On our way to the stable age distribution we computed the intrinsic birth rate. What's the intrinsic death rate?