## Abstract

In this article, we develop a fully integrated and dynamic Bayesian approach to forecast populations by age and sex. The approach embeds the Lee-Carter type models for forecasting the age patterns, with associated measures of uncertainty, of fertility, mortality, immigration, and emigration within a cohort projection model. The methodology may be adapted to handle different data types and sources of information. To illustrate, we analyze time series data for the United Kingdom and forecast the components of population change to the year 2024. We also compare the results obtained from different forecast models for age-specific fertility, mortality, and migration. In doing so, we demonstrate the flexibility and advantages of adopting the Bayesian approach for population forecasting and highlight areas where this work could be extended.

## Introduction

This work is guided by two aims. The first is to have a flexible platform for forecasting populations. Most statistical offices in developed countries utilize data obtained from different sources, including administrative registers, surveys, and censuses with varying levels of quality and measurement. Cohort component models have long been the standard apparatus for producing projections, but wide differences remain in the underlying assumptions, especially regarding the treatment of migration and the level of detail provided. The forecasting approach we have developed is one that can be adapted to include all demographic components of change by age and sex, and provides measures of the accuracy of the forecasts.

As we move away from deterministic population projections to those that provide measures of uncertainty, we believe it is important to integrate the various sources of uncertainty into the modeling framework. The rationale for considering a Bayesian approach is that it offers a natural probabilistic framework to predict future populations. Variability in the data and uncertainties in the parameters and model choice can be explicitly incorporated by using probability distributions, and the predictive distributions follow directly from the probabilistic model applied. The approach also allows the inclusion of expert judgments, together with their uncertainty, in the model framework.

The second aim of this article is to provide a flexible and consistent method for forecasting the age patterns of fertility, mortality, and migration that drive our forecast results. A vast literature focuses on modeling demographic events (e.g., Bongaarts and Bulatao 2000), but approaches for forecasting the age patterns are less abundant. Methods for forecasting migration, in particular, represent a persistent weakness in population forecast models (Bijak 2010).

Our population forecasting model is developed with these two aims in mind. We focus on generalizing and extending the Lee-Carter model for forecasting mortality to age-specific fertility and migration (Lee and Carter 1992), and integrating these into the cohort projection mechanism. One of the contributions of the proposed approach is forecasting age-specific emigration rates and immigration volumes, following the suggestions of Rees (1986:148). Because the age patterns of immigration and emigration are more regular than those observed for net migration, they are better amenable to modeling by using Lee-Carter models.

## Background

### Forecasting the Age Patterns of Fertility, Mortality, and Migration

There is a long history of modeling the age patterns of fertility, mortality, and migration events (Booth 2006). This work has demonstrated the persistent and strong regularities in the age patterns over time and across space (Preston et al. 2001:191–210; Rogers 1986; Rogers and Little 1994) driven by biological and social life course mechanisms (Courgeau 1985). The age regularities exhibited in demographic patterns allow population forecasters to simplify their underlying assumptions and models. Indeed, some forecasts focus on indicator variables, such as the total fertility rate (TFR), life expectancy, or net migration rate, which are then converted into an assumed age distribution (see, e.g., Raftery et al. 2012; Wilson and Bell 2007).

The main approaches for modeling the age patterns of demographic components include the imposition of empirical tables obtained from other countries and historical settings (Coale and Demeny 1966), parametric model schedules (Coale and McNeil 1972; Coale and Trussell 1974; Heligman and Pollard 1980; Knudsen et al. 1993; Rogers 1986; Rogers and Castro 1981; Rogers et al. 1978), relational models (Brass 1974), functional models (De Beer 2011; Hyndman and Booth 2008; Hyndman and Ullah 2007; Lee and Carter 1992), and hierarchical Bayesian models (Czado et al. 2005; Girosi and King 2008). Of these, the most successful and widely used for forecasting future patterns and uncertainty have been the relational and functional models—particularly, the Lee-Carter model for mortality (Booth and Tickle 2008).

Relational and functional models include standards or time-invariant patterns, which are then perturbed based on a set of parameters. Because the shapes of age-specific fertility, mortality, and migration largely remain the same over time, these approaches provide a powerful tool for forecasting. In the Lee-Carter case, single- or five-year age groups are altered on the basis of time series equations. The main critique of the original Lee-Carter approach is that it can produce implausible forecasts for particular age groups (see Girosi and King 2008:38–42). As a consequence, several extensions have been developed to accommodate cohort effects, correlation between sexes, and smoothing (Booth and Tickle 2008; Lee 2000), as well as to forecast fertility rates (Lee 1993).

### Bayesian Population Forecasting

The need to incorporate probabilistic uncertainty into population estimates and forecasts is well known. Probabilistic forecasts have the advantage over variant style projections in that they specify the chances or probability that a particular future population value will be within any given range (Ahlburg and Land 1992; Alho and Spencer 1985, 2005; Bongaarts and Bulatao 2000; Keilman 1990; Lee and Tuljapurkar 1994; Lutz 1996). With variant projections, on the other hand, the user has no idea how likely future population values are, but only that they are plausible scenarios representing the “most likely” and the “extreme” high and low possibilities. However, despite the known advantages of probabilistic forecasts, they have yet to be widely adopted by statistical agencies (Lutz and Goldstein 2004). The reason is that there are many types of uncertainties to consider, and including them in projections is not always straightforward.

According to Alho (1999), probabilistic population forecasting within the Bayesian framework has a tradition dating back over 60 years to the seminal works of Leo Törnqvist and colleagues (Hyppölä et al. 1949). However, it was not until the 1980s that probabilistic methods began entering mainstream demography. These included time series extrapolations, expert-based forecasting, and past error propagation (for a detailed overview of different approaches, see Bijak 2010). Examples of Bayesian models for population forecasting were practically non-existent until the past few years, with the notable exception of Daponte et al. (1997). Recent advances in fast computation and numerical methods have enabled a more widespread use of the Bayesian approach in many fields of application, including population forecasting. At a global level, Raftery et al. (2012) proposed a generic model for all countries of the world, based on aggregate indicators (TFRs and life expectancies) and model life tables. At the opposite end of the data spectrum, Bryant and Graham (2013) suggested a comprehensive Bayesian approach to reconcile different data sources for New Zealand, a country with very good availability of population statistics.

Bryant and Graham’s (2013) approach uses an accounting framework for estimating New Zealand’s current population disaggregated by regions, age, sex and time. It combines various data sources, including vital events registers, censuses, and school and electoral rolls. The model constrains the true values of the demographic components by the population accounting equation. Also, age, time, sex and regional patterns are specified by main effect and two-way interaction terms within a Poisson-gamma model, similar to non-Bayesian approaches of Smith et al. (2010) and Raymer et al. (2011b).

Concurrently, Wheldon et al. (2013) undertook Bayesian estimation and projection of populations to reconstruct past population data. Their approach is based on modeling the three population components—fertility, mortality, and net migration—and accounts for the varying quality of the population figures available from the censuses. Census data are treated as biased estimates of the true unknown population count. Their approach does not provide a systematic modeling of the age profiles and does not account for changing behaviors over time. The information about the model parameters is fed into the model in the form of the informative prior distributions. The component forecasts are inserted in the cohort-component model similar to the one described in the upcoming subsection on the population projection model.

## Methodology

In this section, we first introduce the forecasting model proposed by Lee and Carter (1992) and then describe how it can be extended and applied within a Bayesian framework. The Lee-Carter model was originally designed to forecast age-specific mortality rates with the following specification:
$logμx,t=αx+βxκt+ξx,t,$
(1)
where the logarithm of the age and time-specific mortality rate μx,t is decomposed into an overall age profile, αx, averaged over the entire period under consideration, and age-specific changes in mortality βx. The subscripts x and t denote age and time, respectively. The βx parameter describes how fast the rates decline over time in response to changes in the time-specific effect κt. The error term ξx,t is assumed to be normally distributed with a mean of 0 and a constant variance. To forecast mortality rates into the future, a simple random walk with drift model for κt was proposed:
$κt=ϕ+κt−1+εt.$
(2)
To ensure identifiability of the model parameters, constraints are imposed such that a sum of βx over age is 1 and a sum of κt over time is 0 (Lee and Carter 1992:661). Lee (1993) subsequently proposed a similar model to forecast age-specific fertility for the United States. In that model, several constraints were introduced to represent the prior information on fertility.
In this article, we extend and adapt the Lee-Carter model to create a general framework for forecasting mortality, fertility, emigration and immigration. In this framework, we first assume that count data on a given population component Yx,tg,k follow a Poisson distribution:
$Yx,tg,k~Poissonμx,tg,kRx,tg,k,$
(3)
where g ∈ {D, B, E, I} denotes a component of population change: D represents deaths (mortality), B is births (fertility), E is emigration, I is immigration, sex k = F denotes females, and k = M stands for males. Parameter μx,tg,k denotes a ratio of counts of demographic event scaled to the size of the population exposed to risk of these events, Rx,tg,k. As in the original Lee-Carter model, t and x denote time and age, respectively. For mortality and emigration, the population at risk is the same; for fertility, it consists of women of reproductive age. For immigration, we forecast the counts rather than rates; hence, we assume that Rx,tI,k = 1 for all x, t, and k. Czado et al. (2005) developed the extension of the Lee-Carter model to incorporate Poisson variability of death counts within the Bayesian framework.
Second, we assume that the logarithm of the rate follows a normal distribution:
$logμx,tg,k~Nαxg,k+βxg,kκtg,k+γt−xg,k,τg,k,$
(4)
where αxg,k, βxg,k, and κtg,k represent the same parameters as in the Lee-Carter model, and γt − xg,k denotes a cohort effect. The cohort effect, introduced by Renshaw and Haberman (2006) for mortality, is incorporated in our framework for the sake of generality, but it may be omitted if not required. Throughout this article, N(μ, τ) denotes a normal distribution with a mean of μ and precision (inverse variance) τ. The normal distribution assumed for rates is an extension of the Czado et al. (2005) model. It allows capturing the overdispersion that is not explained by the variability resulting from the Poisson sampling of count data.

Third, for the time-specific effects, κtg,k, we require a time series model, which facilitates the forecasting. This model can be univariate, such as random walk with drift in the original Lee-Carter specification, for each component and each sex; alternatively, the model can be multivariate (e.g., vector autoregression (VAR)), which allows exploring correlations between sexes and components. A time series model is also utilized for a cohort effect, γt − xg,k. In our application, we use a univariate model. However, more general multivariate frameworks can be used.

To ensure identification of the parameters αxg,k, βxg,k, κtg,k, and γt − xg,k, the following constraints are imposed:
$∑x=0zβxg,k=1,κ1g,k=0,γ1g,k=0,$
(5)
where z denotes the oldest age group. These constraints suffice to identify the bilinear model in Eq. (4) as long as there is a clear differentiation in the βx—that is, as long as they are not all equal to 1 / z, in which case the model reduces to the linear age-period-cohort (APC) model. The problem of identification of the period and cohort effects in the APC model has long been discussed in the literature (see Luo 2013, with a comment by Fienberg 2013).

To learn about the model parameters, and specifically about the forecasts of the population components, we use Bayesian inference. Bayes theorem states that the posterior distribution of the model parameters (e.g., forecasts of the age-specific mortality rates) is proportional to the product of the likelihood for the data and the prior distribution. In our approach, Bayesian inference integrates uncertainty from the demographic event (count) data expressed by the Poisson likelihood with weakly informative prior distributions that give preference to the historical data. Subsequently, all four components of population change are combined within the cohort component projection model.

In the following subsections, we present specific adaptations of the preceding framework to forecast mortality, fertility, and migration. We adopt a convention of proposing a very simple model for the data (such as the original Lee-Carter model) with a more general one. Because our extensions of the Lee-Carter model lead to a relatively complex specification of the probabilistic model, the closed forms of the posterior distributions are difficult to obtain analytically. Hence, we sample from the posterior distributions by using the Markov chain Monte Carlo (MCMC) algorithms implemented in the OpenBUGS software (Lunn et al. 2009). The example code used for the simulations for fertility is available in Online Resource 1 (section A.3). The other codes are available from the corresponding author upon request.

### Forecasting Mortality

To forecast the mortality of males and females, we consider the Bayesian version of the original Lee-Carter model, denoted as M1, and a general extension of this model, denoted by M2. The Lee-Carter model M1 is specified as in Eqs. (1) and (2), and the age-specific mortality rates are calculated as μx,tD,k = Yx,tD,k / Rx,tD,k. Model M2 for death counts Yx,tD,k follows Eqs. (3) and (4).

In M1, the time-specific parameters κtk (we drop the superscript D for each parameter for clarity of notation) follow univariate random walk with drift models, as in the original Lee-Carter specification. In M2, κtk for both sexes follow a bivariate vector autoregressive VAR(1) process with drift:
$κtFκtM~MNV2ϕ01ϕ02+ϕ11ϕ12ϕ21ϕ22κt−1Fκt−1M,Tκ,$
(6)
where MVN2 denotes a two-dimensional multivariate normal distribution, with the precision matrix Tκ; ϕij, i, j = 1, 2 are the parameters of the VAR(1) model, with ϕ0j, j = 1, 2 being drift terms. The instantaneous and lagged correlations assumed for males and females reflect the assumption of parallel improvements in health conditions over time.1 Finally, the cohort effect γt − xk for each sex k follows a univariate autoregressive process AR(1) with parameters ψ0k and ψ1k:
$γt−xk~Nψ0k+ψ1kγt−x−1k,τγk.$
(7)

### Forecasting Fertility

For forecasting age-specific fertility rates, we apply a simple version of the Lee (1993) model, here called F1. The extended model, denoted by F2, includes a cohort effect (see also Cheng and Lin 2010; Lee 1993, 2000). The population at risk represents all women of reproductive age.

The time component in F1 follows an ARMA (1,1) process (again we suppress the superscripts B and k because the model relates only to females):
$κt~Nϕ0+ϕ1κt−1−ϕ0+ϕ2ξt−1,τκ,$
(8)
where ξt = κt − ϕ0 − ϕ1t − 1 − ϕ0) − ϕ2ξt − 1.
In Model F2, we use a simple univariate autoregressive process AR(1) for the time component κt and for γt − x:
$κt∼Nϕ0+ϕ1κt−1,τκ,γt−x∼Nψ0+ψ1γt−x−1,τγ.$
(9)

### Forecasting Immigration Counts and Emigration Rates

To forecast immigration counts and emigration rates, we introduce two models: (1) a univariate model, denoted by IE1, which assumes no correlation between emigration and immigration of males and females; and (2) a multivariate model, denoted by IE2, in which we assume correlation between the time parameters κt for both sexes and both directions of migration. In both models, we incorporate smoothing, built in into the prior distributions for the age-specific model parameters αx and βx. Unlike mortality and fertility, there is no clear rationale for including a cohort effect parameter in forecasting immigration and emigration. Also, because immigration does not have an easily defined population at risk, we model the counts, which is a common practice in population projections (McDonald and Kippen 2002; Rees 1986).

In IE1, we assume a random walk without drift for emigration rates and immigration counts for both sexes:
$κtg,k~Nκt−1g,kτκg,k.$
(10)
This simple model leads to forecasts with a constant expectation (as the last observation) and increasing uncertainty.
For IE2, we assume instantaneous correlation in the time parameters for emigration and immigration for both sexes:
$κt~MVN4ϕ0+ϕ1logt+ϕ2t+∘ϕ3κt−1,Tκ,$
(11)
where κt = (κtE, F, κtE, M, κtI, F, κtI, M)′, ϕi = (ϕi1, . . . , ϕi4)′, Tκ is a precision matrix, (∘) denotes element-wise multiplication, and the prime notation (′) denotes transposition. We simplify the model by having the time parameters, κt, depend only on their own lagged values and not on the direction of the flow. This model includes a drift term ϕ0, an autoregressive parameter ϕ3, a logarithmic trend with parameter ϕ1, and a linear trend with parameter ϕ2.

### Prior Distributions

In the absence of any prior information, we suggest using weakly informative distributions that will allow the data to drive the estimation. For the general forecasting model, we propose a set of priors, the hyperparameters of which may differ in particular applications. For both sexes k, the specifications of the prior distributions for the model parameters are as follows:
$αxk~N00.01forallx,β1:z−1k~MVNz−11/z,τβkΨβk,βzk=1−∑i=1z−1βxk,τβk~Γ0.0010.001,ϕij∼N01,foralliandj,ψik∼N01,i=0,1,σk∼U0100,τk=σk−2,σγk∼U0100,τγk=σγk−2,Tκ∼WishartlIl,l,$
(12)
where Γ(a, b) denotes the gamma distribution with a mean of a / b and variance a / b2; U denotes uniform distribution; and Ψβk is a precision matrix for a conditional distribution of βx, given that they sum to 1 (details are presented in Online Resource 1, section A.2). The preceding prior distributions imply weak information a priori about the model parameters. For the age-specific parameters αx and βx, we avoid specifications based on the data (i.e., the empirical Bayes approach) suggested by Czado et al. (2005). For the precision parameters τk and τγk, we follow Gelman’s (2006) suggestions of avoiding conjugate gamma distributions, and we use uniform priors for standard deviations over a large range. The Wishart distribution is a standard prior distribution for precision matrix in the multivariate time series model and is easy to implement in the OpenBUGS software. The subscript l denotes the number of series in the model: males and females for mortality (i.e., l = 2) or male and female migration in both directions (i.e., l = 4). In the univariate model (e.g., for fertility or for a single sex), we replace the Wishart prior with a uniform prior for the standard deviation; that is, σκ ~ U(0, 100), τκ = (σκ)− 2.

For low-quality data, smoothing of the age profile may be required. In this case, we propose a smoothing technique that is embedded in the specification of the prior distributions for parameters αx and βx. Smoothing prevents the artificial age patterns resulting from the sample data from being propagated in the forecasts. Our smoothing method is based on the spatial autoregressive processes (see, e.g., Besag 1986).

Prior densities for smoothing the age-specific parameters αx are constructed in the following way, with the sex and population component superscripts omitted for clarity. For the youngest and oldest age groups, we assume that the mean of the prior distribution depends on the second-youngest and second-oldest group, respectively:
$α0~Nα1,12τα,αz~Nαz−1,12τα.$
(13)
For the remaining age groups, we assume that their means depend on the average of the two neighboring age groups x – 1 and x + 1:
$αx~N12αx−1+12αx+1,τα,x∉0z.$
(14)
The preceding construction of the conditional precisions for each age group ensures that the unconditional precision is constant for all age groups. Finally, we assume a priori that τα = 100, which is based on our assessment of results obtained with various values for this parameter and from visual inspections of the model fits. It implies a moderate degree of smoothing for αx.
For βx, we assume the same pattern of smoothing as for αx, but we derive a distribution conditional on ∑βx = 1. The resulting multivariate normal distribution is
$β1:z−1∼MVNz−11/z,τβΨβ,βz=1−∑x=1z−1βxg,k,τβ~Γ0.000010.00001.$
(15)
This prior distribution is similar to the one in Eq. (12). However, the matrix Ψβ here is derived analogously to Eqs. (13) and (14) by assuming that all elements of βx follow an autoregressive process that in the limit tends to a random walk, but also conditional on ∑βx = 1 (for more details see Online Resource 1, section A.2). The smoothing parameter τβ can be sex-specific and direction-of-flow-specific, or one parameter can be used for all four flows, which allows borrowing of strength. The gamma distribution assumed for this parameter is characterized by a very heavy tail; thus, it allows this parameter to explore regions of large values that lead to “smoother” age profiles. Because the smoothing parameter has a vague prior distribution, the whole smoothing procedure is driven by the data rather than by subjective judgment. However, the results exhibit sensitivity to the specification of the prior for the smoothing parameter. Other prior distributions, such as truncated t or Cauchy, can be used. The degree of smoothing may also be controlled by fixing τβ at some value, which can be found by a grid search, for example.

### Population Projection Model

The results of forecasting the four components of population change—that is, samples from the posterior distributions of mortality, fertility, and emigration rates, as well as immigration counts—are subsequently combined into a cohort component projection model (see Preston et al. 2001:117–137; Rogers 1995). The projection model is specified as
(16)
where Ptk = (P0,tk, …, Pz,tk)′ is vector of midyear population sizes by age and sex k and Itk = (μ0,tI,k, …, μz,tI,k)′ is a vector of forecasted immigration counts. Further, btk = (0, …, b14,tk, …, b45,tk, …, 0) is vector of birth rates,
$stk=s0,tk00⋯00s1,tk0⋯0⋮⋱⋮00⋯sz−2,tk0000⋯0sz−1,tksz,tk$
is a matrix of survivorship rates, and a = 1 / 2.05 is the assumed proportion of female births in the population. Finally, 0 = (0, . . . , 0) is a vector of length z, and O is a matrix of zeros of size (z – 1 × z). The survivorship rates come from the mortality and emigration models (Rogers 1995:104–107):
$sx,tk=1−0.5μx,tD,k+μx,tE,k1+0.5μx+1,tD,k+μx+1,tE,k,forx≠z,$
(17)
$sz,tk=1−0.5μz,tD,k+μz,tE,k1+0.5μz,tD,k+μz,tE,k.$
(18)
They include the standard transformation of the mortality and emigration rates, both of which result in the decrease of the population, into the survival rates. For the last year of age x = z, we assume the standard formula following Rogers (1995:107). Finally, the birth rates are constructed by using the fertility rates obtained from our model for forecasting fertility and survivorship of infants from the mortality forecasting model:
$bx,tk=11+0.5μ0,tD,k12μx,tB,F+sx,tFμx+1,tB,F.$
(19)

### Model Validation and Selection

The models for the population components that underlie the population forecast are selected from the models proposed in the previous sections. The selection process is based on (1) visual evaluation of goodness of fit of the model to the data and the forecasts, (2) ex-post evaluation of the in-sample forecasts of the population components based on the 1975–2000 truncated data set and, where appropriate, (3) the deviance information criterion (DIC) as a formal criterion for model selection (Spiegelhalter et al. 2002).

The DIC is a tool for assessing the goodness-of-fit of a model to the data, which enables selecting the best performing model. It is often considered a generalization of the Akaike information criterion (AIC) for comparing complex Bayesian hierarchical models. It utilizes a deviance of the likelihood evaluated at the mean of the posterior distribution of the likelihood as the goodness-of-fit measure, corrected with the “effective number of parameters” in a model (for definitions and formal derivation, see Spiegelhalter et al. 2002). The requirement for using the DIC is that the posterior distribution is approximately multivariate normal. Selecting the best performing model is similar for the AIC—namely, the lower the value of the criterion, the better fit of the model.

## Illustration: The Case of the United Kingdom

### Data

In this section, we illustrate our forecasting method with the data for the United Kingdom (UK). These data represent a case in which the counts of all population components by single year of age and sex are available but their quality is varying. Because the data on vital events are recorded by the registers, they are considered to be precise and of relatively good quality. Immigration and emigration counts are, however, produced by using the International Passenger Survey (IPS) and include both sampling and nonsampling errors, especially in the age profiles by single year of age.

The data used to produce our forecasts represent the period 1975–2009. The data on mortality rates were obtained from the Human Mortality Database (n.d.). The emigration and immigration counts were obtained from the Office for National Statistics. The data on births were obtained from the Office for National Statistics (England and Wales), Northern Ireland Statistics Research Agency, and National Records of Scotland. The UK midyear population estimate for 2009, used as a baseline for predictions, was also obtained from the Office for National Statistics. Logarithms of single year mortality rates for females and males from 1975 to 2009 are presented in the upper row of Fig. 1. We observe that (1) mortality at all ages, and for both sexes, have been decreasing over time; (2) females have lower mortality than males; and (3) males exhibit considerably higher mortality in the young adult years.

Fertility rates by age of mother are presented in the bottom row of Fig. 1. Over time, we observe a shift from a peak level of fertility at ages 23–26 in 1970 toward one at ages 29–33 years in 2009. The reasons for this shift are related to fertility postponement and a subsequent recuperation. Because of the relatively small counts for very young and very old ages, the data on births were aggregated into age groups under 15 years and 45 years and older. To compute fertility rates, the same female population at risk that was used to calculate the age-specific mortality was applied, except for the age groups under 15 and 45 and older, for which the population at risk was aggregated for ages 12–14 years and 45–50 years, respectively. Further, in our illustration, an implicit assumption was made about fertility—namely, that the rates for boundary ages (i.e., under 15 years and 45 years and older) are applied to the population aged 14 years and 45 years, respectively. However, because these rates are very small, the overall effect is negligible.

The total flows of immigration and emigration from 1975 to 2009 are presented in the top row of Fig. 2. We observe similar trends in male and female migration over time. The immigration levels increased rapidly from the 1990s through around 2005. For emigration, the increase is less noticeable and appears to be more volatile, which may be caused by random sample variation in the underlying data source, the International Passenger Survey. Larger irregularities appear when the data are disaggregated by single year of age, as illustrated for immigration and emigration in the middle and bottom rows, respectively, of Fig. 2 (see also Raymer et al. 2011a).

## Results

In this section, we present the results of forecasting the population components with the models described in the previous section. For each component, we discuss the model’s goodness-of-fit to the data and forecasts of the future patterns, and select the underlying model to be used for the population forecast.

### Forecasts of Mortality

In the first row of Fig. 3, we present the fit of the models M1 and M2 to the 2009 data. It can be observed that the fit of the model M2 with the cohort component reflects the data better than M1, which is especially clear when comparing the fits for life expectancy (third and fourth rows). In particular, M2 is able to reproduce the mortality volatility of the cohort born during the influenza pandemic in 1918–1919. Mortality projected with M1 is lower than that projected with M2 for age groups 0–15 and 35–70, and the age pattern is more uncertain, as depicted in the second row of Fig. 3. For life expectancy (third row), the fit to the data is more uncertain under M1 than M2; in the case of predictions, however, uncertainty is larger under M2. Also, M1 leads to lower predicted life expectancy than M2.

Recent literature has pointed to the importance of the cohort effect in measuring and predicting period mortality rates and the resulting life expectancies (Luy 2010; Luy and Wegner 2009). In particular, cohort effects are likely to stem from the long-lasting effects of early-life events and circumstances on mortality, rather than being a result of whole life trajectories experienced by particular cohorts, as demonstrated in a series of longitudinal studies (e.g., Bengtsson and Mineau 2009; for a general overview and a critical discussion, see also Murphy 2010). An alternative argument for the inclusion of the cohort effect in the model is lifelong processes that might affect mortality, such as smoking (e.g., Doll et al. 2004).

To support our rationale for selecting model M2, we analyze the in-sample forecasts from both models based on the 1975–2000 truncated data set. Model M1 (the original Lee-Carter model) yields forecasts that slightly underpredict the observed life expectancy, which is presented in the bottom plot of Fig. 3. For M2, the posterior distribution is much wider than the results obtained with the full data set. Here, the median predictions seem to be slightly lower than the observed life expectancy, but the observed values are inside the predictive intervals because of larger uncertainty. A comparison of the ex-post forecasts reveals that 59 % of the observed mortality rates for years 2001–2009 fall into the 80 % predictive intervals in the M1 model. For the 95 % predictive interval, 82 % of the observations fall into it. The M2 model performs better; the percentages of the observed mortality rates falling into 80 % and 95 % predictive intervals are 75 % and 89 %, respectively. Because we model mortality rates in M1 and death counts in M2, the DIC cannot be used here to compare both models.

Our life expectancy forecasts can be compared with the official ones prepared by the Office for National Statistics (2011). For 2024, the official predictions of 85.3 for females and 81.6 for males fall inside the 80 % predictive intervals. Median life expectancies are 83.9 and 80.0 under M1, and 85.1 and 80.6 under M2. Hence, the model with cohort effect (M2) leads to slightly lower predictions of life expectancy compared with the official ones, but higher predictions compared with M1.

### Forecasts of Fertility

The age-specific forecasts for fertility are presented in Fig. 4. In the first row, we observe the fit of the models F1 and F2 to the 2009 data. The model with the cohort effect (F2) provides a better fit with lower uncertainty. Also, the 2024 forecast (second row) appears more plausible than the forecast based on the F1 model, which produces an unrealistic median fertility rate of 0.3 for females aged 33–35 years.

The resulting TFRs are presented in the third row of Fig. 4. It is clear from the plots that F2 fits the data better than F1. Moreover, the projected TFR from F1 shows an explosive pattern that we consider unrealistic, with an explosive predicted TFR. Hence, we believe that the pattern of gradual diminishing of the recently increasing TFR produced by F2 better reflects our expectations about future fertility in the UK.

The in-sample forecasts of the fertility rates confirm our rationale for choosing F2 as the foundation of the population forecast. Again, F2 appears to fit the data better (see the fourth row of Fig. 4). The resulting forecasts of TFR under F2 seem to be more uncertain than those of F1. However, F1 misses the decline in early 2000s. These results are confirmed by the ex-post analysis of the fertility rates. For F1, 58 % of observed fertility rates fall into the 80 % predictive interval, and 75% fall into the 95 % predictive intervals; for F2, the percentage of data falling into respective predictive intervals are 62 % and 74 %. The official forecast of the TFR used by the Office for National Statistics (2011) is 1.84, and it falls inside the 80 % predictive interval of our 2014 forecast under F2. Our median TFR forecast for 2024 is 1.12.

The DIC cannot be used to compare F1 and F2 because different types of data are used in the models. Nevertheless, the ex-post analysis of the in-sample forecasts, as well as the visual assessment of the results, clearly point to the model with the cohort effect included. This rationale is supported by the vast demographic literature on the quantum and tempo effects in fertility (Bongaarts and Feeney 1998). In particular, we refer to the recent postponement and subsequent recuperation of fertility in many developed countries, where the cohort effects are the most profound (see, e.g., Sobotka et al. 2011). In our results, slightly declining but still uncertain fertility rates may indicate a possibility of yet another period of postponement.

#### Forecasts of Emigration and Immigration

To forecast emigration rates and immigration counts for the UK we need to reflect on the quality of the survey data on migration and their implications for modeling. To overcome the irregularities described in the earlier Data section, we modified the models described in the section Forecasting Immigration Counts and Emigration Rates to account for rounding to the nearest thousand and included smoothing in the model. Detailed equations are presented in Online Resource 1 (section A.1).

The results of forecasting emigration rates and immigration counts for females are summarized in Figs. 5 and 6, respectively (the similar patterns for males are not shown because of space constraints). In the first row, we present the IE1 and IE2 forecasts for 2009. We observe that both models fit the data reasonably well. As expected, the univariate random walk model in IE1 leads to stable forecasts over time, with increasing uncertainty for both emigration rates and immigration counts. The drift term and log-linear trend incorporated in the IE2 lead to ever-increasing immigration and emigration.

The DIC leads to choosing the IE2 model over IE1 in the case of the full sample predictions: the DIC is 25,484 for IE1 and 25,480 for IE2. In the case of the truncated data set, however, the DIC prefers the IE1 with random walk (18,928 vs. 18,930). This is supported by the visual inspection of the in-sample predictions of mean emigration rates and total immigration counts (fourth rows of Figs. 5 and 6, respectively). This result is not surprising because the patterns in the migration data changed substantially after the year 2000. Hence, different models may be more suitable for both data sets. In terms of predictive coverage, there is almost no difference between IE1 and IE2. Under IE1, 37 % of observed migration fall into the 80 % predictive interval, and 48 % fall into the 95 % predictive interval; for IE2, the percentages of data falling into respective predictive intervals are 37 % and 47 %. The rather small percentages of the observed values falling into the predictive intervals ought not to be surprising due to the irregularities observed in the data.

### Population Forecasts

A final step in forecasting population is combining the population components within the cohort component projection model. As an illustration, we select models M2 for mortality, F2 for fertility, and IE2 for migration.

The age composition of the predicted UK population in 2024 is presented in the first row of Fig. 7. Forecasts of the total population for females and males are presented in the second row. We observe that the age profile of the 2024 population is shaped mostly by future migration and, to a lesser extent, fertility. The largest uncertainty concerns the youngest population, as well as the population aged 20–45, for both males and females. These findings are in line with Keyfitz’s (1981) observation on the plausible limits of population forecasting, which were set to about 20 years ahead. We also expect that the number of the elderly persons will be larger in 2024 but the population aged 20–45 will be most numerous.

Our forecast can be compared with the official deterministic projections for 2024 prepared by the Office for National Statistics (2011). The dashed line in Fig. 7 represents the principal projection, whereas the dotted lines are low and high population scenarios. The main drivers of differences between these projections are assumptions about fertility and migration. We observe that the population aged 20–35 is substantially smaller in all scenarios of the Office for National Statistics projections compared with our forecast. This results from the assumption that net migration stays constant at the levels observed in recent years: at the level of 200,000, 140,000, and 260,000 persons annually in the principal, low, and high scenarios, respectively (Office for National Statistics 2011). Moreover, our results suggest that the uncertainty about the age profile of future population is considerably larger than it is reflected in the deterministic projections based on scenarios.

As far as the total population size is concerned, the forecast indicates that there will be only 49,000 fewer females than males in 2024, whereas the difference was more than 1 million in 2009 and was 1.5 million in favor of females in 1975. This is most likely due to the larger proportion of male migration and a gradual closing of the life expectancy gap between the sexes. The median size of the 2024 population is 70.8 million, which is around 9 million larger than the population size observed in 2009. In the principal projection for 2024 prepared by the Office for National Statistics (2011), the predicted total population is 69.0 million (i.e., 1.8 million lower), whereas the low and high scenarios are 66.4 and 71.1 million, respectively. However, only high scenario falls into our 80 % predictive interval for the total population. Also, the population predicted by the Office for National Statistics seems to be growing more slowly than in our forecast. This results from the aforementioned more-conservative prediction of international net migration.

As a further validation step in population forecasting, we present an in-sample population forecast in Fig. 8, based on the data truncated in year 2000. In general, the data can be truncated at various points to make the validation procedure more robust. The predictive intervals for the age profiles in 2009 match the reported figures reasonably well. Differences are observed for early childhood ages, for which the model underpredicts the fertility increases in the first decade of the 2000s, and for the young adult ages, caused by the underprediction of migration (especially for females). The EU enlargement in May 2004 resulted in a faster increase than anticipated by the historical data for both emigration and immigration levels, which explains a large part of the underprediction.

## Discussion

In this article, we have demonstrated that extension of the Lee-Carter model can serve as a general platform for estimating age schedules of the four demographic components of population. We then combined these components into a single forecast by means of a cohort-component projection model. We also explored the correlation of each of the components in time, as well as between sexes and components (for emigration and immigration), which is embedded in our extension of the Lee-Carter method. For emigration and immigration, we provided a tool for smoothing irregularities in the data. This tool, however, can be easily extended to fertility or mortality. Finally, we have illustrated the use of the forecasting model on the UK’s data.

This research makes two contributions to the literature. The first is the development of a new approach for integrating demographic components to provide stochastic population forecasts by age and sex. The Bayesian approach that we adopted accounts for the uncertainties embedded in births, deaths, and emigration and immigration, as well as across age and sexes. We show that the same general framework of the Lee-Carter approach for modeling age and sex patterns of mortality and fertility can be coherently applied to model corresponding patterns of migration. Irregularities in the data, such as those observed for the UK, can also be accounted for within the model.

The second contribution is the application of the approach to a situation of relatively good yet imperfect data availability. In this way, we position our work between the generic global approach with far fewer data requirements, which has been proposed by Raftery et al. (2012), and a specific multiple-data situation discussed by Bryant and Graham (2013). Where possible, population forecasting should follow a bottom-up approach, in which the age-specific rates of the demographic components are utilized. The rates describe the underlying processes more comprehensively than summary aggregates, such as TFRs or life expectancies, in the top-down approach.

Further research should explore other models for forecasting age patterns of demographic components, such as the functional models developed by Hyndman and Booth (2008). Analogously, various specifications for the time component models (such as ARIMA or VAR models of higher order) should be investigated. Next, the underlying models of components for the population forecast can be selected by using various techniques, of which Bayesian model averaging (Raftery et al. 1997) seems to be most appealing. In this way, the model uncertainty would be accounted for coherently. Finally, the uncertainty of the baseline population size used for projections could be incorporated into the projection. We believe this work provides a strong foundation for such extensions.

## Acknowledgments

We gratefully acknowledge a grant from the Economic and Social Research Council Centre for Population Change (Grant RES-625-28-0001). We thank the Editor and four anonymous reviewers for their valuable comments and suggestions on this work; Jack Baker, Heather Booth, Jennifer M. Ortman, Adrian Raftery and participants of the annual meeting of the Population Association of America 2013, New Orleans.

## Note

1

Li and Lee (2005) proposed an alternative approach to include correlations between sexes; they added a commonality factor to Eq. (1).

## References

Ahlburg, D. A., & Land, K. C. (Eds.). (
1992
). Population forecasting: Introduction.
International Journal of Forecasting
,
8
,
289
299
.
Alho, J. M. (
1999
,
August
).
On probabilistic forecasts of population and their uses
. Paper presented at the 52nd Session of the International Statistical Institute, Helsinki, Finland.
Alho, J. M., & Spencer, B. D. (
1985
).
Uncertain population forecasting
.
Journal of the American Statistical Association
,
80
,
306
314
. 10.1080/01621459.1985.10478113
Alho, J. M., & Spencer, B. D. (
2005
).
Statistical demography and forecasting
.
New York, NY
:
Springer
.
Bengtsson, T., & Mineau, G. P. (
2009
).
Early-life effects on socio-economic performance and mortality in later life: A full life-course approach using contemporary and historical sources
.
Social Science & Medicine
,
68
,
1561
1564
. 10.1016/j.socscimed.2009.02.012
Besag, J. E. (
1986
).
On the statistical analysis of dirty pictures
.
Journal of the Royal Statistical Society: Series B
,
48
,
259
302
.
Bijak, J. (
2010
).
Forecasting international migration in Europe: A Bayesian view
.
Dordrecht, The Netherlands
:
Springer
.
Bongaarts, J., & Bulatao, R. A. (Eds.). (
2000
).
Beyond six billion: Forecasting the world’s population
.
Washington, DC
:
.
Bongaarts, J., & Feeney, G. (
1998
).
When is a tempo effect a tempo distortion?
.
Genus
,
66
,
1
15
.
Booth, H. (
2006
).
Demographic forecasting: 1980 to 2005 in review
.
International Journal of Forecasting
,
22
,
547
581
. 10.1016/j.ijforecast.2006.04.001
Booth, H., & Tickle, L. (
2008
).
Mortality modelling and forecasting: A review of methods
.
Annals of Actuarial Science
,
3
,
3
43
. 10.1017/S1748499500000440
Brass, W. (
1974
).
Perspectives in population prediction: Illustrated by the statistics of England and Wales
.
Journal of the Royal Statistical Society: Series A
,
137
,
532
570
. 10.2307/2344713
Bryant, J. R., & Graham, P. J. (
2013
).
Bayesian demographic accounts: Subnational population estimation using multiple data sources
.
Bayesian Analysis
,
8
(
2
),
1
32
.
Cheng, P. C. R., & Lin, E. S. (
2010
).
Completing incomplete cohort fertility schedules
.
Demographic Research
,
23
(
article 9
),
223
256
. 10.4054/DemRes.2010.23.9
Coale, A., & Demeny, P. (
1966
).
Regional model life tables and stable populations
.
Princeton, NJ
:
Princeton University Press
.
Coale, A., & McNeil, D. R. (
1972
).
The distribution by age of the frequency of first marriage in a female cohort
.
Journal of the American Statistical Association
,
67
,
743
749
. 10.1080/01621459.1972.10481287
Coale, A., & Trussell, J. (
1974
).
Model fertility schedules: Variations in the age structure of childbearing in human populations
.
Population Index
,
40
,
185
258
. 10.2307/2733910
Courgeau, D. (
1985
).
Interaction between spatial mobility, family and career life-cycle: A French survey
.
European Sociological Review
,
1
,
139
162
.
Czado, C., Delwarde, A., & Denuit, M. (
2005
).
Bayesian Poisson log-bilinear mortality projections
.
Insurance: Mathematics and Economics
,
36
,
260
284
.
Daponte, B., Kadane, J., & Wolfson, L. (
1997
).
Bayesian demography: Projecting the Iraqi Kurdish population, 1977–1990
.
Journal of the American Statistical Association
,
92
,
1256
1267
.
De Beer, J. (
2011
).
A new relational method for smoothing and projecting age-specific fertility rates: TOPALS
.
Demographic Research
,
24
(
article 18
),
409
454
. 10.4054/DemRes.2011.24.18
Doll, R., Peto, R., Boreham, J., & Sutherland, I. (
2004
).
Mortality in relation to smoking: 50 years’ observations on male British doctors
.
British Medical Journal
,
328
,
1519
. 10.1136/bmj.38142.554479.AE
Fienberg, S. E. (
2013
).
Cohort analysis’ unholy quest: A discussion
.
Demography
,
50
,
1981
1984
. 10.1007/s13524-013-0251-z
Gelman, A. (
2006
).
Prior distributions for variance parameters in hierarchical models
.
Bayesian Analysis
,
1
,
515
533
. 10.1214/06-BA117A
Girosi, F., & King, G. (
2008
).
Demographic forecasting
.
Princeton, NJ
:
Princeton University Press
.
Heligman, L., & Pollard, J. H. (
1980
).
The age pattern of mortality
.
Journal of the Institute of Actuaries
,
107
,
49
80
. 10.1017/S0020268100040257
Human Mortality Database
. (n.d.).
University of California, Berkeley (USA), and Max Planck Institute for Demographic Research (Germany)
Hyndman, R. J., & Booth, H. (
2008
).
Stochastic population forecasts using functional data models for mortality, fertility and migration
.
International Journal of Forecasting
,
24
,
323
342
. 10.1016/j.ijforecast.2008.02.009
Hyndman, R. J., & Ullah, M. S. (
2007
).
Robust forecasting of mortality and fertility rates: A functional data approach
.
Computational Statistics and Data Analysis
,
51
,
4942
4956
. 10.1016/j.csda.2006.07.028
Hyppölä, J., Tunkelo, A., & Törnqvist, L. (
1949
).
Suomen väestöä, sen uusiutumista ja tulevaa kehitystä koskevia laskelmia
[Calculations concerning the population of Finland, its renewal and future development] (Tilastollisia tiedonantoja, 38).
Helsinki
:
Statistics Finland
.
Keilman, N. (
1990
).
Uncertainty in national population forecasting: Issues, backgrounds, analyses, recommendations
.
Amsterdam, The Netherlands
:
Swets & Zeitlinger
.
Keyfitz, N. (
1981
).
The limits of population forecasting
.
Population and Development Review
,
7
,
579
593
. 10.2307/1972799
Knudsen, C., McNown, R., & Rogers, A. (
1993
).
Forecasting fertility: An application of time series methods to parameterized model schedules
.
Social Science Research
,
22
,
1
23
. 10.1006/ssre.1993.1001
Lee, R. D. (
1993
).
Modeling and forecasting the time series of US fertility: Age distribution, range, and ultimate level
.
International Journal of Forecasting
,
9
,
187
202
. 10.1016/0169-2070(93)90004-7
Lee, R. D. (
2000
).
The Lee-Carter method for forecasting mortality, with various extensions and applications
.
North American Actuarial Journal
,
4
,
80
93
. 10.1080/10920277.2000.10595882
Lee, R. D., & Carter, L. R. (
1992
).
Modelling and forecasting U.S. mortality
.
Journal of the American Statistical Association
,
87
,
659
671
.
Lee, R. D., & Tuljapurkar, S. (
1994
).
Stochastic population projections for the United States: Beyond high, medium and low
.
Journal of the American Statistical Association
,
89
,
1175
1189
. 10.1080/01621459.1994.10476857
Li, N., & Lee, R. (
2005
).
Coherent mortality forecasts for a group of populations: An extension of the Lee-Carter method
.
Demography
,
42
,
575
594
. 10.1353/dem.2005.0021
Lunn, D., Spiegelhalter, D., Thomas, A., & Best, N. (
2009
).
The BUGS project: Evolution, critique, and future directions
.
Statistics in Medicine
,
28
,
3049
3067
. 10.1002/sim.3680
Luo, L. (
2013
).
Assessing validity and application scope of the intrinsic estimator approach to the age-period-cohort problem
.
Demography
,
50
,
1945
1967
. 10.1007/s13524-013-0243-z
Lutz, W. (Ed.). (
1996
).
The future population of the world: What can we assume today?
.
London, UK
:
Earthscan
.
Lutz, W., & Goldstein, J. R. (
2004
).
Introduction: How to deal with uncertainty in population forecasting?
.
International Statistical Review
,
72
,
1
4
. 10.1111/j.1751-5823.2004.tb00219.x
Luy, M. (
2010
).
Tempo effects and their relevance in demographic analysis
.
Comparative Population Studies
,
35
,
415
446
.
Luy, M., & Wegner, C. H. (
2009
).
Conventional versus tempo-adjusted life expectancy: Which is the more appropriate measure for period mortality?
.
Genus
,
65
,
1
28
.
McDonald, P., & Kippen, R. (
2002
).
Projecting future migration levels: Should rates or numbers be used?
.
People and Place
,
10
,
82
83
.
Murphy, M. (
2010
).
Reexamining the dominance of birth cohort effects on mortality
.
Population and Development Review
,
36
,
365
390
. 10.1111/j.1728-4457.2010.00334.x
Office for National Statistics
. (
2011
).
National Population Projections, 2010-based reference volume: Series PP2
.
Titchfield, UK
:
Population Projections Unit, Office for National Statistics
.
Preston, S. H., Heuveline, P., & Guillot, M. (
2001
).
Demography: Measuring and modeling population processes
.
Oxford, UK
:
Blackwell
.
Raftery, A. E., Li, N., Ševcíková, H., Gerland, P., & Heilig, G. K. (
2012
).
Bayesian probabilistic population projections for all countries
.
Proceedings of the National Academy of Sciences of the USA
,
109
,
13915
13921
. 10.1073/pnas.1211452109
Raftery, A. E., Madigan, D., & Hoeting, J. A. (
1997
).
Bayesian model averaging for linear regression models
.
Journal of the American Statistical Association
,
92
,
179
191
. 10.1080/01621459.1997.10473615
Raymer, J., Abel, G. J., Disney, G., & Wiśniowski, A. (
2011a
).
Improving estimates of migration flows to Eurostat
(CPC Working Paper 15).
Southampton, UK
:
ESRC Centre for Population Change, University of Southampton
.
Raymer, J., De Beer, J., & Van der Erf, R. (
2011
).
Putting the pieces of the puzzle together: Age and sex-specific estimates of migration amongst countries in the EU/EFTA, 2002–2007
.
European Journal of Population
,
27
,
185
215
. 10.1007/s10680-011-9230-5
Rees, P. H. (
1986
).
Choices in the construction of regional population projections
. In R. I. Woods & P. H. Rees (Eds.),
Population structures and models: Developments in spatial demography
(pp.
126
159
).
London, UK
:
George Allen and Unwin
.
Renshaw, A. E., & Haberman, S. (
2006
).
A cohort-based extension to the Lee-Carter model for mortality reduction factors
.
Insurance: Mathematics and Economics
,
58
,
556
570
.
Rogers, A. (
1986
).
Parameterized multistate population dynamics and projections
.
Journal of the American Statistical Association
,
81
,
48
61
. 10.1080/01621459.1986.10478237
Rogers, A. (
1995
).
Multiregional demography: Principles, methods and extensions
.
Chichester, UK
:
John Wiley & Sons
.
Rogers, A., & Castro, L. J. (
1981
).
Model migration schedules
(Research Report 81-30).
Laxenburg, Austria
:
International Institute for Applied Systems Analysis
.
Rogers, A., & Little, J. S. (
1994
).
Parameterizing age patterns of demographic rates with the multiexponential model schedule
.
Mathematical Population Studies
,
4
,
175
195
. 10.1080/08898489409525372
Rogers, A., Raquillet, R., & Castro, L. J. (
1978
).
Model migration schedules and their applications
.
Environment and Planning A
,
10
,
475
502
. 10.1068/a100475
Smith, P. W. F., Raymer, J., & Giulietti, C. (
2010
).
Combining available migration data in England to study economic activity flows over time
.
Journal of the Royal Statistical Society: Series A
,
173
,
733
753
. 10.1111/j.1467-985X.2009.00630.x
Sobotka, T., Zeman, K., Lesthaeghe, R., Frejka, T., & Neels, K. (
2011
).
Postponement and recuperation in cohort fertility: Austria, Germany and Switzerland in a European context
.
Comparative Population Studies
,
36
,
417
452
.
Spiegelhalter, D. J., Best, N. G., Carlin, B. P., & Van der Linde, A. (
2002
).
Bayesian measures of model complexity and fit
.
Journal of the Royal Statistical Society: Series B
,
64
,
585
639
.
Wheldon, M. C., Raftery, A. E., Clark, S. J., & Gerland, P. (
2013
).
Reconstructing past populations with uncertainty from fragmentary data
.
Journal of the American Statistical Association
,
108
,
96
110
. 10.1080/01621459.2012.737729
Wilson, T., & Bell, M. (
2007
).
Probabilistic regional population forecasts: The example of Queensland, Australia
.
Geographical Analysis
,
39
,
1
25
. 10.1111/j.1538-4632.2006.00693.x