## 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

_{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:

_{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.

*Y*

_{x,t}

^{g,k}follow a Poisson distribution:

*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,t}

^{g,k}denotes a ratio of counts of demographic event scaled to the size of the population exposed to risk of these events,

*R*

_{x,t}

^{g,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

*R*

_{x,t}

^{I,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.

_{x}

^{g,k}, β

_{x}

^{g,k}, and κ

_{t}

^{g,k}represent the same parameters as in the Lee-Carter model, and

*γ*

_{t − x}

^{g,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, κ_{t}^{g,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 − x}^{g,k}. In our application, we use a univariate model. However, more general multivariate frameworks can be used.

_{x}

^{g,k}, β

_{x}

^{g,k}, κ

_{t}

^{g,k}, and

*γ*

_{t − x}

^{g,k}, the following constraints are imposed:

*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,t}^{D,k} = *Y*_{x,t}^{D,k} / *R*_{x,t}^{D,k}. Model M2 for death counts *Y*_{x,t}^{D,k} follows Eqs. (3) and (4).

_{t}

^{k}(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, κ

_{t}

^{k}for both sexes follow a bivariate vector autoregressive VAR(1) process with drift:

*MVN*

_{2}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 − x}

^{k}for each sex

*k*follows a univariate autoregressive process AR(1) with parameters ψ

_{0}

^{k}and ψ

_{1}

^{k}:

### 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.

*B*and

*k*because the model relates only to females):

*ξ*

_{t}= κ

_{t}− ϕ

_{0}− ϕ

_{1}(κ

_{t − 1}− ϕ

_{0}) − ϕ

_{2}

*ξ*

_{t − 1}.

*κ*

_{t}and for

*γ*

_{t − x}:

### 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).

**κ**

_{t}= (κ

_{t}

^{E, F}, κ

_{t}

^{E, M}, κ

_{t}

^{I, F}, κ

_{t}

^{I, 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

*k*, the specifications of the prior distributions for the model parameters are as follows:

*a*,

*b*) denotes the gamma distribution with a mean of

*a*/

*b*and variance

*a*/

*b*

^{2};

*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).

_{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:

*x*– 1 and

*x*+ 1:

_{α}= 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}.

_{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

**Ψ**

_{β}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

**P**

_{t}

^{k}= (

*P*

_{0,t}

^{k}, …,

*P*

_{z,t}

^{k})′ is vector of midyear population sizes by age and sex

*k*and

**I**

_{t}

^{k}= (μ

_{0,t}

^{I,k}, …, μ

_{z,t}

^{I,k})′ is a vector of forecasted immigration counts. Further,

**b**

_{t}

^{k}= (0, …,

*b*

_{14,t}

^{k}, …,

*b*

_{45,t}

^{k}, …, 0) is vector of birth rates,

*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):

*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:

### 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

**Open Access**This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.