This computing log illustrates the use of the cohort component method in population projections.
We will use the data from Sweden used in Box 6.1 (page 125) in the
textbook, available as file sweden93.dat in the datasets
section to save you some typing. We read the population counts, the life
table, and the fertility rates. The last age group is 85+ and
L for that age is actually time lived after age 85. We also
express person-years per birth, dividing by l0.
. infile age p L f using /// > https://grodri.github.io/datasets/sweden93.dat, clear (18 observations read) . mata ───────────────────────────────────────────────── mata (type end to exit) ────── : p93 = st_data(., "p") : L = st_data(., "L")/100000 : f = st_data(., "f") : sum(p93) 4397428 : end ────────────────────────────────────────────────────────────────────────────────
> library(dplyr)
> sw <- read.table("https://grodri.github.io/datasets/sweden93.dat",
+ header=FALSE)
> names(sw) <- c("age", "p93", "L", "f")
> sw <- mutate(sw, L = L/100000)
> summarize(sw, sum(p93))
sum(p93)
1 4397428
I have written a function to compute a Leslie matrix given two vectors representing person-years lived 5Lx (with the convention that the radix is one and the last value is Tx) and the maternity rates 5mx. We source and then list the code.
. capture mata mata drop Leslie()
. quietly do leslie.mata
. type leslie.mata
// Leslie matrix for population projection
// G. Rodríguez / 28feb2017
//
mata:
real matrix Leslie(real vector L, real vector m) {
n = length(L)
M = J(n,n,0)
// lower diagonal has survivorship ratios
for (i=1; i < n; i++) {
M[i+1,i] = L[i+1]/L[i]
}
M[n,n-1] = M[n,n] = L[n]/(L[n-1]+L[n])
// first row has net maternity contributions
for(i=1; i < n; i++) {
if(m[i]==0 & m[i+1]==0) continue
M[1,i] = L[1]*(m[i] + m[i+1]*L[i+1]/L[i])/2
}
if (m[n] > 0) M[1,n] = L[1]*m[n]
return(M)
}
end
> source("leslie.R")
> Leslie
function (L, m)
{
n = length(L)
M = matrix(0, n, n)
for (i in 1:(n - 1)) {
M[i + 1, i] <- L[i + 1]/L[i]
}
M[n, n - 1] <- M[n, n] <- L[n]/(L[n - 1] + L[n])
for (i in 1:(n - 1)) {
if (m[i] != 0 | m[i + 1] != 0) {
M[1, i] <- L[1] * (m[i] + m[i + 1] * L[i + 1]/L[i])/2
}
}
if (m[n] > 0)
M[1, n] <- L[1] * m[n]
M
}
The function computes the survivorship ratios and stores them in the
lower diagonal. For example M[2,1] = L[2]/L[1] is the
probability of surviving from age 0-4 to 5-9, or
5L5/5L0 in standard
demographic notation.
The only tricky bit is the last (open) age interval, where we combine
the last two groups and project then using
L[n]/(L[n-1] + L[n]). In Sweden we combine ages 80-84 and
85+ and use T85/T80 =
T85/(5L80+T80) as the
survival ratio. (As noted in class, if we had T90 we
would use a slightly different procedure. The textbook describes that
procedure on page 121, but uses the combined projection in Box 6.1.)
The rest of the calculation computes the average fertility rate for
each age group and then survives the resulting births to age 0-4. For
example women 15-19 are exposed to the rates at 15-19 and 20-24, with
the latter discounted by the probability of surviving to 20-24, so we
average m[4] and m[5]*L[5]/L[4].
The rate would be multiplied by 5, the width of the period, and by
the probability of surviving from birth to age 0-4, which is
L[1]/5. The 5’s cancel, so we don’t include them. (As noted
in class, there are two ways of computing births; here we focus on the
women and average the rates, an approach best suited for computing the
Leslie matrix. The textbook also describes focusing on the rates and
averaging the numbers of women exposed to each rate.)
For Sweden we have fertility rates, so we divide by 2.05 to obtain maternity rates (female births) before calling our function. Once we have a Leslie matrix, projection is easy
. mata:
───────────────────────────────────────────────── mata (type end to exit) ──────
: M = Leslie(L, f/2.05)
: p98 = M * p93
: p03 = M * p98
: round( (p93, p98, p03) )
1 2 3
┌────────────────────────────┐
1 │ 293395 293574 280121 │
2 │ 248369 293189 293368 │
3 │ 240012 248251 293049 │
4 │ 261346 239833 248066 │
5 │ 285209 261015 239529 │
6 │ 314388 284787 260629 │
7 │ 281290 313782 284238 │
8 │ 286923 280463 312859 │
9 │ 304108 285576 279147 │
10 │ 324946 301731 283344 │
11 │ 247613 320974 298043 │
12 │ 211351 243039 315045 │
13 │ 215140 205109 235861 │
14 │ 221764 204944 195388 │
15 │ 223506 204793 189260 │
16 │ 183654 194419 178141 │
17 │ 141990 142324 150666 │
18 │ 112424 131768 141960 │
└────────────────────────────┘
: sum(p98), sum(p03)
1 2
┌─────────────────────────────┐
1 │ 4449569.727 4478712.273 │
└─────────────────────────────┘
: end
────────────────────────────────────────────────────────────────────────────────
> M <- Leslie(sw$L, sw$f/2.05)
> sw <- mutate(sw,
+ p98 = M %*% p93,
+ p03 = M %*% p98)
> round(select(sw, age, p93, p98, p03), 0)
age p93 p98 p03
1 0 293395 293574 280121
2 5 248369 293189 293368
3 10 240012 248251 293049
4 15 261346 239833 248066
5 20 285209 261015 239529
6 25 314388 284787 260629
7 30 281290 313782 284238
8 35 286923 280463 312859
9 40 304108 285576 279147
10 45 324946 301731 283344
11 50 247613 320974 298043
12 55 211351 243039 315045
13 60 215140 205109 235861
14 65 221764 204944 195388
15 70 223506 204793 189260
16 75 183654 194419 178141
17 80 141990 142324 150666
18 85 112424 131768 141960
> select(sw, p93, p98, p03) %>% colSums()
p93 p98 p03
4397428 4449570 4478712
These results agree exactly with Box 6.1 (part 2) in the text.
While we have the Leslie matrix handy, we can compute the intrinsic rate of growth and the stable age distribution by simply obtaining the first eigenvalue and eigenvector of the projection matrix.
This is easy to do in Mata, as long as we define the output vector and matrix before. Here I use empty matrices.
. mata
───────────────────────────────────────────────── mata (type end to exit) ──────
: values = J(0,0,.)
: vectors = J(0,0,.)
: eigensystem(M, vectors, values)
: values[1]
1.00111253
: log(values[1])/5
.000222383
: stable = Re(vectors[,1]/sum(vectors[,1]))
: stable
1
┌───────────────┐
1 │ .0619116748 │
2 │ .0617994881 │
3 │ .0617013816 │
4 │ .0615869203 │
5 │ .061440548 │
6 │ .0612814302 │
7 │ .0610952699 │
8 │ .0608479527 │
9 │ .0604950315 │
10 │ .0599554177 │
11 │ .0591567528 │
12 │ .0579994637 │
13 │ .056223923 │
14 │ .0534997799 │
15 │ .0493506342 │
16 │ .0428804537 │
17 │ .0331935317 │
18 │ .035580346 │
└───────────────┘
: end
────────────────────────────────────────────────────────────────────────────────
> e <- eigen(M) > e$values[1] [1] 1.001113+0i > log(e$values[1])/5 [1] 0.0002223831+0i > arv <- abs(Re(e$vectors[,1])) > stable <- arv/sum(arv) > stable [1] 0.06191167 0.06179949 0.06170138 0.06158692 0.06144055 0.06128143 [7] 0.06109527 0.06084795 0.06049503 0.05995542 0.05915675 0.05799946 [13] 0.05622392 0.05349978 0.04935063 0.04288045 0.03319353 0.03558035
The first eigenvalue is 1.0011, so if the 1993 fertility and
mortality rates were to stay constant the population would eventually
grow 0.11 percent every 5 years, or 0.022 percent per year. The first
eigenvector, divided by its sum, gives us the proportionate age
distribution. In general the eigen values and eigen vectors can be
complex numbers, and both Stata and R store them as such, although for a
Leslie matrix they are always real. The function Re() takes
the real part of a complex number.
We return to stable populations in a separate handout, but there is a lot you could do with the Swedish data. For one thing you might try projecting the population for 100 years of so (20 periods) to verify that it becomes stable. You may also try plotting the current and stable equivalent age distributions to see what that says about population momentum.