Probabilistic Projection of Subnational Total Fertility Rates

TL;DR

This paper develops a probabilistic method for projecting subnational total fertility rates that aligns with UN standards and demonstrates superior predictive performance across diverse countries.

Contribution

It introduces a region-specific scaling method for TFR projections that improves accuracy and supports the Watkins hypothesis of within-country convergence over time.

Findings

The scaling method outperforms other approaches in predictive accuracy.

It effectively captures within-country TFR correlation.

Supports the Watkins hypothesis over the last 50 years.

Abstract

We consider the problem of probabilistic projection of the total fertility rate (TFR) for subnational regions. We seek a method that is consistent with the UN's recently adopted Bayesian method for probabilistic TFR projections for all countries, and works well for all countries. We assess various possible methods using subnational TFR data for 47 countries. We find that the method that performs best in terms of out-of-sample predictive performance and also in terms of reproducing the within-country correlation in TFR is a method that scales the national trajectory by a region-specific scale factor that is allowed to vary slowly over time. This supports the hypothesis of Watkins (1990, 1991) that within-country TFR converges over time in response to country-specific factors, and extends the Watkins hypothesis to the last 50 years and to a much wider range of countries around the world.

Tables

Table 1. Table 1: Out-of-sample validation of probabilistic subnational TFR projections over 1995-2010. MAE is mean absolute error. CRPS is continuous ranked probability score. The 80 % percent 80 80\% and 95 % percent 95 95\% columns refer to the percentage of the observations that fell within their prediction interval. The marginal TFR was validated on 3199 values; the average TFR was validated on 137 values. The Scale-AR(1) parameters were ϕ = 0.898 italic-ϕ 0.898 \phi=0.898 and σ = 0.0533 𝜎 0.0533 \sigma=0.0533 .

	Marginal TFR					Average TFR
	MAE	bias	CRPS	80%	95%	MAE	bias	CRPS	80%	95%
Scale-AR(1)	$0.205$	$-0.088$	$-0.147$	$82.0$	$96.3$	$0.172$	$-0.117$	$-0.127$	$82.5$	$95.6$
1d-BHM	$0.228$	$-0.067$	$-0.167$	$75.1$	$90.1$	$0.169$	$-0.101$	$-0.123$	$71.5$	$89.1$
1d-BHM (indep)	$0.228$	$-0.071$	$-0.167$	$75.2$	$89.8$	$0.169$	$-0.103$	$-0.142$	$38.7$	$50.4$
Scale	$0.220$	$-0.106$	$-0.156$	$76.2$	$92.2$	$0.182$	$-0.136$	$-0.133$	$78.8$	$95.6$
Persistence	$0.365$	$-0.305$	$-0.365$	–	–	$0.334$	$-0.303$	$-0.334$	–	–

Table 2. Table 2: Out-of-sample validation of TFR projections over 1990-2010. The marginal TFR was validated on 4144 values; the average TFR was validated on 180 values. The Scale-AR(1) parameters were ϕ = 0.910 italic-ϕ 0.910 \phi=0.910 and σ = 0.0513 𝜎 0.0513 \sigma=0.0513 .

	Marginal TFR					Average TFR
	MAE	bias	CRPS	80%	95%	MAE	bias	CRPS	80%	95%
Scale-AR(1)	$0.323$	$-0.209$	$-0.234$	$70.0$	$84.8$	$0.278$	$-0.192$	$-0.202$	$73.3$	$87.2$
1d-BHM	$0.344$	$-0.207$	$-0.260$	$64.5$	$79.1$	$0.284$	$-0.187$	$-0.214$	$60.0$	$76.1$
1d-BHM (indep)	$0.344$	$-0.209$	$-0.260$	$64.6$	$79.6$	$0.284$	$-0.190$	$-0.242$	$26.7$	$41.7$
Scale	$0.333$	$-0.230$	$-0.245$	$65.6$	$82.0$	$0.291$	$-0.215$	$-0.210$	$72.2$	$87.2$
Persistence	$0.590$	$-0.538$	$-0.590$	–	–	$0.519$	$-0.476$	$-0.522$	–	–

Table 3. Table 3: Out-of-sample validation of average TFR projections over 1995-2010, validated on 137 values.

	MAE	bias	CRPS	80%	95%
1d-BHM (indep.)	$0.169$	$-0.103$	$-0.142$	$38.7$	$50.4$
1d-BHM (1.)	$0.170$	$-0.101$	$-0.124$	$69.3$	$89.1$
1d-BHM (2.)	$0.167$	$-0.102$	$-0.122$	$74.5$	$89.8$
1d-BHM (3.)	$0.168$	$-0.104$	$-0.127$	$62.0$	$80.3$
1d-BHM (4.)	$0.167$	$-0.105$	$-0.127$	$62.8$	$81.0$
1d-BHM (5.)	$0.167$	$-0.106$	$-0.124$	$69.3$	$84.7$
1d-BHM (6.)	$0.167$	$-0.106$	$-0.125$	$68.6$	$83.9$
1d-BHM (7.)	$0.168$	$-0.104$	$-0.125$	$65.7$	$83.2$
1d-BHM (8.)	$0.169$	$-0.101$	$-0.123$	$71.5$	$89.1$
1d-BHM (9.)	$0.169$	$-0.100$	$-0.123$	$70.1$	$89.1$
1d-BHM (10.)	$0.168$	$-0.102$	$-0.123$	$70.8$	$89.1$
1d-BHM (11.)	$0.167$	$-0.103$	$-0.122$	$71.5$	$89.8$
Scale-AR(1)	$0.172$	$-0.117$	$-0.127$	$82.5$	$95.6$
Scale	$0.182$	$-0.136$	$-0.133$	$78.8$	$95.6$
Persistence	$0.334$	$-0.303$	$-0.334$	–	–

Table 4. Table 4: Out-of-sample validation of average TFR projections over 1990-2010, validated on 180 values.

	MAE	bias	CRPS	80%	95%
1d-BHM (indep.)	$0.284$	$-0.190$	$-0.242$	$26.7$	$41.7$
1d-BHM (1.)	$0.284$	$-0.187$	$-0.213$	$60.0$	$77.2$
1d-BHM (2.)	$0.279$	$-0.190$	$-0.210$	$63.3$	$78.3$
1d-BHM (3.)	$0.284$	$-0.192$	$-0.217$	$55.6$	$71.1$
1d-BHM (4.)	$0.283$	$-0.192$	$-0.216$	$55.6$	$70.0$
1d-BHM (5.)	$0.284$	$-0.194$	$-0.215$	$57.8$	$76.1$
1d-BHM (6.)	$0.284$	$-0.194$	$-0.216$	$57.8$	$74.4$
1d-BHM (7.)	$0.284$	$-0.191$	$-0.214$	$56.7$	$72.2$
1d-BHM (8.)	$0.284$	$-0.187$	$-0.214$	$60.0$	$76.1$
1d-BHM (9.)	$0.284$	$-0.187$	$-0.213$	$60.0$	$76.1$
1d-BHM (10.)	$0.285$	$-0.188$	$-0.215$	$60.0$	$76.7$
1d-BHM (11.)	$0.284$	$-0.189$	$-0.214$	$61.1$	$77.2$
Scale-AR(1)	$0.278$	$-0.192$	$-0.202$	$73.3$	$87.2$
Scale	$0.291$	$-0.215$	$-0.210$	$72.2$	$87.2$
Persistence	$0.519$	$-0.476$	$-0.522$	–	–

Equations

$$f_{c,t+1}-{\mu }_c={\rho }_c(f_{c,t}-{\mu }_c)+{\epsilon }_{c,t},\text{with}{\epsilon }_{c,t}\stackrel{iid}{\sim }N(0,{\sigma }_{\epsilon }^2).$$

$$f_{r_c,t,i}={\alpha }_{r_c}{f}_{c,t,i},$$

$${\alpha }_{r_c}=f_{r_c,t=P}/f_{c,t=P}.$$

$${\alpha }_{r_c,t}-1=\varphi ({\alpha }_{r_c,t-1}-1)+{\epsilon }_{r_c,t},\text{with}{\epsilon }_{r_c,t}\stackrel{iid}{\sim }N(0,{\sigma }_c^2).$$

$$\varphi$$

$${\sigma }_c^2$$

$$\sigma$$

$${\mathbf{\epsilon }}_{\mathbf{t}}\sim N(0,{\Sigma }_t={\mathbf{\sigma }}_{\mathbf{t}}^{\prime }\mathbf{A}{\mathbf{\sigma }}_{\mathbf{t}}),$$

$$\mathbf{A}=\frac{\bar{T}-1}{\bar{T}}\tilde{\mathbf{A}}+\frac{1}{2\bar{T}}$$

$${\alpha }_{r_c,t}-1=\varphi ({\alpha }_{r_c,t-1}-1)+{\epsilon }_{r_c,t},\text{with}{\epsilon }_{r_c,t}\stackrel{iid}{\sim }N(0,{\sigma }_c^2).$$

$$\underset{t\to \infty }{\lim }\text{Var}({\alpha }_{r_c,t})\le {\text{Var}}_{q\in R_c}({\alpha }_{q,t=P}),$$

$${\sigma }_c^2=\min \{{\sigma }^2,(1-{\varphi }^2){\text{Var}}_{r\in R_c}({\alpha }_{r,t=P})\}$$

$$\Delta {\alpha }_{r_c,t}={\alpha }_{r_c,t}-{\alpha }_{r_c,t-1}.$$

$$\underset{t\to \infty }{\lim }\text{Var}({\alpha }_{r_c,t})$$

$$\underset{t\to \infty }{\lim }\text{Var}(\Delta {\alpha }_{r_c,t})$$

$$\varphi$$

$${\sigma }^2$$

$$\text{Var}({\alpha }_{r_c,t})$$

$$\text{Var}(\Delta {\alpha }_{r_c,t})$$

$$\text{SD}({\alpha }_{r_c,t})$$

$$\text{SD}(\Delta {\alpha }_{r_c,t})$$

$$\varphi$$

$$\sigma$$

$${\mathbf{\epsilon }}_{\mathbf{t}}\sim N(0,{\Sigma }_t={\mathbf{\sigma }}_{\mathbf{t}}^{\prime }\mathbf{A}{\mathbf{\sigma }}_{\mathbf{t}}).$$

$${\tilde{A}}_{r,s}=\frac{{\sum }_T{e}_{r,t}{e}_{s,t}}{\sqrt{{\sum }_T{e}_{r,t}^2}\cdot \sqrt{{\sum }_T{e}_{s,t}^2}}.$$

$$B$$

$$C$$

$$\hat{\mathbf{A}}=\frac{\bar{T}-1}{\bar{T}}{{\tilde{\mathbf{A}}}^{\ast }}_{[\ge 0]}+\frac{1}{2\bar{T}},$$

$$\hat{\mathbf{A}}=\frac{T-1}{T}\mathbf{B}+\frac{1}{2T}\mathbf{J},$$

$$x^{\prime }\hat{\mathbf{A}}x=\frac{T-1}{T}{x}^{\prime }\mathbf{B}x+\frac{1}{2T}{x}^{\prime }\mathbf{J}x.$$

$$x^{\prime }\mathbf{J}x=(\sum _{i=1}^n{x}_i{)}^2\ge 0.$$

$${\tilde{A}}_{r,s}=\frac{1/T{\sum }_T{e}_{r,t}{e}_{s,t}}{\sqrt{1/T{\sum }_T{e}_{r,t}^2}\cdot \sqrt{1/T{\sum }_T{e}_{s,t}^2}}$$

$${\hat{A}}_{r,s}=\frac{{\int }_0^1\rho {\left(\frac{1}{\sqrt{1-{\rho }^2}}\right)}^T\exp \left[-\frac{1}{2(1-{\rho }^2)}\left[S{S}_r-2\rho S{S}_{r,s}+S{S}_s\right]\right]d\rho }{{\int }_0^1{\left(\frac{1}{\sqrt{1-{\rho }^2}}\right)}^T\exp \left[-\frac{1}{2(1-{\rho }^2)}\left[S{S}_r-2\rho S{S}_{r,s}+S{S}_s\right]\right]d\rho },$$

$$\text{MAE}=\frac{1}{RT}\sum _{c=1}^C\sum _{r=1}^{R_c}\sum _{t=1}^T|f_{r_c,t}-{\hat{f}}_{r_c,t}|.$$

Full text

Keywords: Total fertility rate, Subnational projections, Autoregressive model, Bayesian hierarchical model, Scaling model, Correlation.

1. 1 Introduction

2. 2 Data

3. 3 Review of the national Bayesian hierarchical model

4. 4 Methods for subnational projections

5. 4.1 Scale Method

6. 4.2 Scale-AR(1)

7. 4.3 One-directional BHM

8. 4.3.1 Correlation between regions

9. 5 Results

10. 5.1 TFR projections

11. 5.2 Out-of-sample predictive validation

12. 6 Discussion

13. A Appendix: Methods

14. A.1 Estimation of the Scale-AR(1) parameters

15. A.2 Estimating the Error Correlations

16. A.3 Out of Sample Validation Measures

1 Introduction

The United Nations Population Division issued official probabilistic population projections for all countries for the first time in 2015 (United Nations,, 2015), using the methodology described by Raftery et al., (2012). One of the key components of the projection methodology is a Bayesian hierarchical model for the total fertility rate (TFR) in all countries (Alkema et al.,, 2011; Raftery et al.,, 2014; Fosdick and Raftery,, 2014).

Population projections for subnational administrative units, such as provinces, states, counties, regions or départements (hereafter all referred to simply as regions), are of great interest to national and local governments for planning, policy and decision-making (Rayer et al.,, 2009). A common current practice is to generate subnational projections deterministically by scaling national projections (U.S. Census Bureau,, 2016). Specifically, the US Census Bureau provides a workbook for users to generate subnational TFR projections for up to 32 regions. The method requires the user to enter an ultimate TFR level (lower asymptote) to which the regional TFR converges, and a deterministic projection of the national TFR. The subnational TFR is then projected in such a way that it approaches the target TFR with the same rate as the national TFR approaches this target. This method does not yield probabilistic projections.

In this paper we try to address one aspect of the problem, namely probabilistic subnational projections of TFR. Methods for probabilistic subnational projections have been developed for individual countries or parts of countries (Smith and Sincich,, 1988; Tayman et al.,, 1998; Rees and Turton,, 1998; Gullickson and Moen,, 2001; Gullickson,, 2001; Lee et al.,, 2003; Smith and Tayman,, 2004; Wilson and Bell,, 2007; Rayer et al.,, 2009; Raymer et al.,, 2012; Wilson,, 2013); for a review see Tayman, (2011). Our ultimate goal is to extend the UN method for probabilistic projections for all countries to a method for subnational probabilistic projections that is consistent across countries and works well for all regions of all countries.

We contrast two broad approaches to subnational probabilistic projection of TFR. One approach is a direct extension of the UN method (Alkema et al.,, 2011) to subnational data, effectively treating the country in the same way the UN model treats the world, and treating the regions in the same way the UN model treats the countries. Borges, (2015) proposed an approach along these lines for the provinces of Brazil.

The other approach is motivated by the observation of Watkins (1990, 1991) that within-country variation in TFR in Europe decreased over the period of the fertility transition, between 1870 and 1960. This observation has been confirmed for a more recent period for the German-speaking countries (Basten et al.,, 2012), to some extent for India (Arokiasmy and Goli,, 2012; Wilson et al.,, 2012), while the evidence is more equivocal for the United States (O’Connell,, 1981). Watkins posited that this was due to increased integration of national markets, expansion of the role of the state, and nation-building in the form of linguistic standardization over this period. Calhoun, (1993) argued that, of these three mechanisms, only linguistic standardization clearly supported her argument. However, some support for the importance of the role of the nation state for fertility is provided by the fact that nation states have specific and different policies aimed at affecting fertility rates (Tomlinson,, 1985; Chamie,, 1994), and some of these policies have been shown to be effective (Kalwij,, 2010; Luci-Greulich and Thévenon,, 2013). Note that Klüsener et al., (2013) investigated subnational convergence of non-marital fertility in Europe in recent decades, and found that within-country variation increased, in contrast with the trends noted by other authors for overall fertility. Here we consider only overall fertility.

One question is then whether the direct extension of the UN method for countries to the subnational context adequately accounts for this tendency of TFR to converge within countries over time. Note that this extension of the UN method does predict within-country convergence of fertility rates over time during the fertility transition; the question is whether it adequately accounts for this convergence.

To investigate this question, we consider a different general approach, which starts from the national probabilistic projections produced by the UN method, and then scales them for each region by a scaling factor that varies stochastically, but stays relatively constant. This induces more within-country correlation than the direct extension of the UN method. It could be viewed as a probabilistic extension of the method currently used by the U.S. Census Bureau. It is also related to the method of Wilson, (2013), but with some significant differences.

We apply these methods to subnational data on total fertility for 47 countries over the period 1950–2010. We compare our two approaches and several variants in terms of out-of-sample predictive performance. The results shed some light on the Watkins hypothesis of increasing within-country correlation, as well providing some guidance on how to carry out subnational probabilistic TFR projection.

Note that there is a substantial literature on convergence of fertility rates in different countries to one another, with different conclusions argued for (Wilson,, 2001, 2004; Reher,, 2004, 2007; Dorius,, 2008; Wilson,, 2011). Our work here has implications for within-country fertility convergence, but is agnostic about fertility convergence between countries, and so does not have implications for global fertility convergence, for example.

The paper is organized as follows. We first describe the data used in this study and review the model for national probabilistic projections. We then introduce our proposed methodology for subnational probabilistic projections, and present the results. The paper concludes with a discussion.

2 Data

We use subnational data on the TFR for 47 countries (13 in the Americas, nine in the Asia-Pacific region, and 25 in Europe), corresponding to 1,092 regions for the period 1950–2010, collected by the United Nations Population Division. Each country analyzed had a population over one million and a national average TFR below 2.5 in 2010–2015. The geographical level selected for each country was the one with available data for the longest comparable time series. The dataset covers 4.9 billion people. Fig. 1 shows the numbers of regions for each country, which range from two for Slovenia to 96 for France. The data include countries from all the inhabited continents except Africa. The data sources are shown in Appendix Table LABEL:tab:datasource.

Fig. 2 shows an example of the data for four countries (USA, India, Brazil and Sweden). It illustrates that the data vary with respect to the correlation between regions. It also shows that the data started later than 1950 for some regions. In the figure, the national TFR from United Nations, (2013) is shown as a black curve.

3 Review of the national Bayesian hierarchical model

Our starting point for developing a methodology for subnational projections is the probabilistic model for projecting national TFR proposed by Alkema et al., (2011), which has now been adopted by the UN for its official projections. We start by summarizing the main ideas of this Bayesian hierarchical model (BHM). More detail can be found in Alkema et al., (2011) and Raftery et al., (2014).

The model is based on standard fertility transition theory (e.g. Hirschman 1994) , and is compatible with almost all versions of this in the literature. It distinguishes three phases in the evolution of a country’s fertility over time, depicted in the left panel of Fig. 3 for the example of Denmark. Phase I (grey dots) precedes the beginning of the fertility transition and is characterized by high fertility that is stable or increasing. This phase is not modeled as all or nearly all countries have completed this phase. During Phase II, or the transition phase (red dots in the figure), fertility declines from high levels to below the replacement level of 2.1 children per woman. Phase III is the post-fertility transition period (blue dots), during which fertility fluctuates at low levels, possibly recovering towards the replacement level.

To model the fertility declines in each five-year period during Phase II, a double logistic decline function is used. An example of this function is shown in the right panel of Fig. 3. The function is parametrized by a set of country-specific parameters that define the shape of the country’s decline curve. Those parameters are drawn from a world distribution. The resulting BHM is estimated using Markov chain Monte Carlo (MCMC).

Phase III is modeled using a Bayesian hierarchichal first-order autoregressive, or AR(1), process of the form:

[TABLE]

It implies that fertility for country $c$ has a country-specific long-term mean, $\mu_{c}$, and autoregressive parameter, $\rho_{c}$, which are assumed to be drawn from a world distribution. The parameters of this world distribution in turn have a joint prior distribution, thus defining a three-level hierarchical model, where the three levels are the observation, the country and the world. The resulting model is again estimated by MCMC.

The process of estimating Phase II and Phase III parameters results in a set of country-specific decline curves and a set of country-specific AR(1) parameter pairs. Unlike decline curves which can be estimated for all countries, not all countries have experienced Phase III, in which cases the country-specific long-term means and autoregressive parameters cannot be estimated. In such cases, the “world” means and autoregressive parameters are used. The estimated parameters are then used to generate a set of future TFR trajectories yielding probabilistic TFR projections for all countries of the world.

4 Methods for subnational projections

Ideally, we seek a method for generating probabilistic subnational TFR projections that reflects the literature and theory of fertility transitions, is based on the national methodology used by the UN and described above, works well for all countries, is as simple if possible, and yields correlations between regions that are similar to the correlations in the observed data.

We first describe a simple Scale method that provides an initial probabilistic extension of methods used by the U.S. Census Bureau and other national agencies. This simple approach works well from many points of view, but it does not allow for the possibility of crossovers between regions, whereas in fact these do happen. We therefore elaborate this model to allow the scale factor to change stochastically, but slowly over time, yielding the so-called Scale-AR(1) method. Finally we describe a quite different approach, called the one-directional BHM, which directly generalizes the national approach to the subnational context, allowing regions to vary more freely within a country.

4.1 Scale Method

We start with a simple intuitive scale method where, for each trajectory from the probabilistic projection, the regional TFR is simply a product of the simulated national TFR and a time-independent but region-specific scale factor.

Let $f_{c,t,i}$ denote the national TFR projection for country $c$ at time $t$ from trajectory $i$, simulated from its posterior distribution as described above. We model $f_{r_{c},t,i}$, the TFR for region $r_{c}$ of country $c$ at time $t$ in the $i$-th trajectory, by

[TABLE]

where $\alpha_{r_{c}}$ denotes the regional scaling factor derived from the last observed (present) time period denoted by $P$:

[TABLE]

Note that $\alpha_{r_{c}}$ is the same for all trajectories. This method yields a set of regional trajectories $f_{r_{c},t,i}$ and thus yields probabilistic projections of the regional TFRs, $f_{r_{c},t}$.

Our numerical experiments, described below, indicated that this simple method performed surprisingly well. However, it also has a serious drawback. Scaling by a constant factor yields a perfect correlation, i.e. it does not allow for the possibility of crossovers between regions over time. However, such crossovers do happen, and the scale method says that they are impossible, which is not fully satisfactory.

4.2 Scale-AR(1)

To avoid this drawback, and modify the scale method so as to allow for the possibility of crossovers, we propose a variation of the simple Scale method where we model the regional scale factor using a first-order autoregressive, or AR(1), process:

[TABLE]

The regional TFR $f_{r_{c},t,i}$ is then derived as in (1) with the additional lower bound restriction, $f_{r_{c},t,i}>0.5$.

This model implies that the scaling factor will fluctuate around one in the long term. Regardless of its initial value, it will converge to a distribution that is centered around one, and the rate of convergence is determined by the $\phi$ parameter. We use the following settings for the model parameters, estimated from the data for all 47 countries available:

[TABLE]

where $P$ again denotes the present time period and $R_{c}$ denotes the set of regions of country $c$. The minimum restriction in (5) ensures that the variation of $\alpha_{\cdot,t}$ across regions is not larger than the variation in the last observed time period, in line with the Watkins hypothesis and the long-term observed data. Details of how these parameters were estimated are given in the appendix.

This method is related to the method proposed by Wilson, (2013), but there are some significant differences that are discussed in the Discussion section.

4.3 One-directional BHM

Next, we consider an extension of the world three-level BHM which is depicted in Fig. 4. The three levels of this model are the world level, country level and time point or observation level. (The time point level is not shown in the figure). In the world version, information from all countries is combined into the world level, which in turn influences the country level, yielding a two-directional BHM. The prior distribution of the hyperparameters is vague for most parameters, but also reflects expert knowledge in some cases. The model yields a posterior distribution of the world parameters, and the country-specific parameters, which is then used to generate the national projections.

Our extension has a similar setup, but moves down by one level of geography and works in one direction only. Thus the top level of our national model is the country, the next level is the region, and the bottom level is the time point. The upper level of our model corresponds to the country level of the world model, that is, we carry over the country-specific posterior from a world simulation and use it as the distribution of the hyperparameters in our national model (red arrow in Fig. 4). On the lower level, data from all regions of a country are handled individually. The estimation of the regional parameters is informed by the hyperparameters, but the regional level does not influence the country level of the model. The resulting regional posterior distribution is used to project subnational TFR.

Note that many countries do not have historical data on Phase III because they have not yet reached this stage, and so in these cases the country posterior is the same as the world posterior. As a result, all regions of those countries inherit the “world” Phase III parameters.

4.3.1 Correlation between regions

For aggregating TFR over sets of regions, for example for deriving country’s averages, it is important to capture correlation in model errors between regions, as was done by Fosdick and Raftery, (2014) for capturing correlation between countries.

We will model the forecast errors as follows:

[TABLE]

where $\boldsymbol{\sigma_{t}}$ is a vector consisting of the forecast standard deviations for each region. For Phase II this is the standard deviation of the errors in the double logistic model and for Phase III it is the standard deviation of the AR(1) model. In (7), $\boldsymbol{A}$ is a matrix where each element $A_{r,s}$ corresponds to the correlation between the model errors of region $r$ and $s$ over all time periods.

Let $f_{r,t}$ denote the observed TFR for region $r$ at time $t$. We denote by $e_{r,t}$ the normalized forecast error, namely the forecast error divided by its standard deviation. The normalized forecast error $e_{r,t}$ is estimated as follows:

• •

Phase II: For each value $g_{r,t,i}$ of a double logistic (DL) trajectory $i$ and the standard deviation of DL $\sigma_{r,i}$ take $d_{r,t,i}=(f_{r,t}-g_{r,t,i})/\sigma_{ri}$. Then $e_{r,t}$ is the mean of $d_{r,t,i}$ over $i$.

• •

Phase III: For each value $h_{r,t,i}$ of a phase III trajectory $i$ and the standard deviation of these trajectories $\sigma_{\varepsilon,r,i}$, take the difference $d_{r,t,i}=(f_{r,t}-h_{r,t,i})/\sigma_{\varepsilon,r,i}$. Then, $e_{r,t}$ is the mean of $d_{r,t,i}$ over $i$.

We define the correlation matrix $\boldsymbol{A}$ as

[TABLE]

where $\boldsymbol{\tilde{A}}$ is a truncated correlation matrix made positive definite, and $\bar{T}$ is the average number of time periods per region. Here $\boldsymbol{A}$ has an approximate Bayesian interpretation as an approximation of the posterior mean with a uniform $U[0,1]$ on the correlations. Note that $\boldsymbol{A}$ is positive definite. The appendix contains details of the method, as well as other methods for deriving $\boldsymbol{A}$ that we have experimented with.

5 Results

We now compare results from the three methods described in the previous section. All three methods depend on a national BHM simulation. We used a simulation that was used to produce the official UN TFR projections in WPP2012. Our version has 2,000 TFR trajectories for each country and was produced using the bayesTFR R package (Ševčíková et al.,, 2011).

For the Scale-AR(1) method, for each region $r_{c}$ we set the initial scaling factor to $\alpha_{r_{c},P}=f_{r_{c},P}/f_{c,P}$ with $P$ being the last observed time period. Then we produced projections of $\alpha_{r_{c},t}$ for $t>P$ using (3). Finally, (1) was applied, as in the case of the simple Scale method, using each of the 2,000 TFR trajectories for country $c$ as $f_{c,t,i}$. This yielded 2,000 regional TFR trajectories.

For the one-directional BHM (1d-BHM), we ran the regional BHM while using the country posterior from the national BHM simulation. Then we projected 2,000 regional TFR trajectories using a sample of the regional posterior parameters. We explored two versions of this model, one that accounts for correlation between regions’ error terms and one that does not, the latter denoted by “1d-BHM (indep)”.

5.1 TFR projections

We are interested in the marginal predictive distribution of future TFR for each region. We are also interested in how reasonable the joint distributions of the trajectories between regions are. Fig. 5 shows one randomly selected trajectory for all regions of Sweden for various methods. In the top panel the Scale-AR(1) method was used. It can be seen that all trajectories closely follow the corresponding national trajectory (black dashed line), while allowing for occasional crossovers. This creates a similar pattern to that seen in the observed data (to the left from the dotted vertical line). The simple Scale method (not shown in the figure) yield trajectories perfectly parallel to the national trajectory with no crossovers.

The bottom two panels of Fig. 5 show results from the one-directional BHM method. In the middle panel we accounted for correlation between regions, whereas in the bottom panel the regions’ error terms were considered independent. As can be seen, this method does not yield trajectories that closely parallel the national one. Furthermore, if correlation is not taken into account, there are many more crossovers between regions than are typically seen in the past data.

All the 47 countries in our dataset show the same pattern in terms of the differences between the methods. In Fig. 6 we selected three countries for which one trajectory obtained via the Scale-AR(1) method is shown for each region (as in the top panel of Fig. 5). As in the case of Sweden, the trajectories are highly correlated and closely follow the national trajectory.

In Fig. 7 we show the predictive median and 80% prediction interval (red) for three regions of India from the Scale-AR(1) method (the corresponding national projection is shown in grey). They represent three different types of regions found across all countries. The first type (in the left panel, Assam, India) is a region with a current TFR that is very close to the national TFR. In such a case, the regional projection mostly overlaps with the national projection, with a slightly larger prediction interval. The black dotted line in the figure shows the median projection resulting from the simple Scale method. This would also be very close to the national median for regions of this type.

Uttar Pradesh in the center is a type of region where current TFR is substantially higher than the national TFR. The underlying AR(1) process causes the median projection of such a region to converge to the national median in the long term, thus decreasing the gap between them. If simple scaling were applied, that gap would remain constant, resulting in much higher projections of TFR for the region.

Finally, Goa on the right, with its current TFR well below the national one, is projected to increase on average, again yielding a smaller gap between the national and regional medians. Here simple scaling results in much lower projections.

5.2 Out-of-sample predictive validation

We validated our methodology via predictive out-of-sample experiments, one for predicting the period 1995–2010, and one for predicting the period 1990-2010. We first assessed the various methods in terms of average predictive performance over all regions of the 47 countries. To assess their performance for predicting aggregates (and hence, for example, in capturing the between-region correlations), we further assess the predictions of the average TFR across the regions of each country.

For both time periods considered, we removed the data points that corresponded to the time period to be predicted, reestimated the models, generated probabilistic projections with the various methods, and compared the projections with the observed data points. The results are shown in Table 1 for 1995-2010 and Table 2 for 1990-2010. The measures in the left part of each table (Marginal TFR) were derived by comparing the probabilistic projections of TFR for all regions to their observed values. The quantities in the right part of the tables (Average TFR) are derived by comparing a TFR averaged over all regions of each country with the observed average TFR for each country.

For comparison purposes, we also added the Persistence method, in which the TFR stays at the same level over time and so the forecast for all future time periods is equal to the last observed value. While this could be viewed as a straw man forecast, persistence forecasts have been found to perform surprisingly well in many forecasting contexts, and so it is worth making this comparison.

In the tables, the mean absolute error (MAE), the bias and the continuous ranked probability score (CRPS) (Hersbach,, 2000; Gneiting and Raftery,, 2007) are reported. The coverages of the 80% and 95% intervals are also reported. The coverage of a prediction interval is defined as the proportion of the time that the truth lies in the interval. We wish the coverage to be close to the nominal level. Thus, for example, ideally the coverage of the 80% interval would be close to 80%.

The appendix gives details of the derivation of these metrics. For MAE and bias, the smaller the absolute value the better. For the two coverage columns an ideal method would match the numbers to the corresponding percentage. The CRPS is a combination of an error-based and a variation-based measure, and thus we give it a high weight when selecting the best method. In this case, a better method corresponds to a larger value of CRPS.

For the marginal TFR, the Scale-AR(1) method performs best in terms of CRPS, MAE and coverage. The simple Scale method comes in second. However, we would not recommend using the simple Scale method because it produces trajectories that are unrealistic in that they do not allow the possibility of crossovers between regions, as mentioned previously. Note that by design, the Scale-AR(1) method yields larger uncertainty than the simple Scale method, which in this case translates to a better coverage and CRPS. The Scale method includes only the uncertainty from the national BHM model, whereas the Scale-AR(1) method has in addition the uncertainty included in the AR(1) process. There is essentially no difference between the 1d-BHM with and without correlation for the marginal TFR. This is expected, as the correlation plays a role only in aggregated indicators.

For the average TFR, the Scale-AR(1) and 1d-BHM have similar performance in terms of CRPS (one is better in Table 1, the other in Table 2). However, Scale-AR(1) has consistently better coverage. Here we see a big difference in coverage between the two versions of 1d-BHM, which does not have good performance if correlation between regions is not taken into account. The good performance of the Scale-AR(1) method suggests that it is accounting adequately for between-region spatial correlation.

6 Discussion

We have developed several methods for subnational probabilistic projection of TFR, and applied them to data from 47 very diverse countries. All the methods take the national projections from the UN method as their starting point. We found that all the methods we propose performed well in terms of out-of-sample predictive performance, and outperformed a simple baseline persistence method.

In the best method, the national trajectories are scaled by a region-specific scaling factor which itself is allowed to vary stochastically but slowly over time. One competing method treats the regions in the same way as countries are treated in the UN’s BHM, but this does not yield enough within-country correlation. Even when we introduce between-region correlation into this model, it still does not have enough within-country correlation overall.

We have compared several different methods, but there are still others in the literature. Rayer et al., (2009) considered ex-post assessment of predictive uncertainty for U.S. counties, extending the national ex-post approach of Keyfitz, (1981) and Stoto, (1983) to the subnational context. Raymer et al., (2012) used a vector autoregressive model for crude birth rates in three regions of England. While these methods may work well for developed countries that have had low fertility for an extended period, they do not capture the systematic variation in fertility decline rates among higher-fertility countries documented by Alkema et al., (2011), and so they may not be so appropriate for our goal here, of developing a method applicable to countries at all levels of the fertility transition.

The extant method closest to our preferred Scale-AR(1) method is one proposed by Wilson, (2013), who also proposed scaling a national TFR forecast by a region-specific scale factor that varies according to an AR(1) model, and applied it to Sidney, Australia. However, there are several differences between the Scale-AR(1) method we propose here, and Wilson’s approach for TFR. The national TFR forecast used by Wilson is based on an AR(1) process centered around an externally-specified main forecast. As discussed, this may not carry over well to higher-fertility countries. Our method, in contrast, is centered around the probabilistic forecast from the UN’s BHM, which is designed to work well for countries at all fertility levels and includes uncertainty about national projections. Also, in our method the model is statistically estimated, while in Wilson’s approach the parameters are adjusted manually.

Our preferred Scale-AR(1) method does not incorporate spatially-indexed between-region correlation. Instead, spatial correlation is modeled by a strong country effect. Our 1-d BHM method did incorporate spatial correlation in the variant that includes between-region correlation estimated from the data (especially methods 8–11 described in the Appendix section on estimating the error correlations). However, this did not allow us to include enough between-region correlation. This may be because within-country correlation seems to be dominated by a strong country effect rather than spatially indexed correlation, as can be seen for example for Sweden in Fig. 5. This is also shown by the good calibration of the Scale-AR(1). Thus we feel it is likely that adding additional spatial correlation would not substantially improve fit of the model to the data at hand.

In addition to providing guidance for subnational projections, our results give insight into how subnational fertility evolves in a modern context. They suggest that there is substantial within-country correlation and convergence. This confirms the observations and hypotheses of Watkins (1990, 1991) for Europe to 1960. It further extends them from just Europe to a range of countries from around the world, and indicates that, broadly speaking, similar patterns continue to hold a half-century later.

Acknowledgements:

This work was supported by NICHD grants R01 HD054511 and R01 HD070936.

Appendix A Appendix: Methods

A.1 Estimation of the Scale-AR(1) parameters

Here we give details on estimating parameters of the Scale-AR(1) model.

The model is based on an AR(1) process for region-specific scale factors $\alpha_{r_{c},t}$ centered at one, namely

[TABLE]

We impose the restriction that the scale factors not diverge indefinitely over time. We implement this by requiring that $\sigma^{2}_{c}$ is such that

[TABLE]

where $P$ denotes the present time period and $R_{c}$ denotes the set of regions in country $c$. This yields

[TABLE]

We are interested in estimating the country- and region-independent parameters $\phi$ and $\sigma$. We know from the observed data that the standard deviation of $\alpha_{r_{c},t}$ declines as TFR declines, which is also in line with the theoretical expectations of Watkins (1990, 1991).

Thus we need to find asymptotic values for those parameters.

Let $\Delta\alpha_{r_{c},t}$ denote the first order differences over time, namely

[TABLE]

Then

[TABLE]

Equations (13) and (14) imply that

[TABLE]

Assuming a normal distribution of $\alpha_{r_{c},t}$ we can write

[TABLE]

From the observed data we know that both $|\alpha_{r_{c},t}-1|$ and $|\Delta\alpha_{r_{c},t}|$ decline as TFR declines (see Fig. 8). The nonparametrically estimated conditional expectation of $|\alpha_{r_{c},t}-1|$ given TFR reaches a minimum, as a function of TFR, of $0.09475$ at TFR $=1.768$, as shown by the dotted lines in Fig. 8. At this level of TFR, the nonparametrically estimated value of $E(|\Delta\alpha_{r_{c},t}|)$ is 0.03678, which is close to its minimum. Taking the mean of $\alpha_{r_{c},t}$ to be 1, and using the fact that the standard deviation of a normal random variable is $\sqrt{\pi/2}$ times its mean absolute deviation, we find that

[TABLE]

Substituting the values from Equations (19) and (20) for Var$(\alpha_{r_{c},t})$ and Var$(\Delta\alpha_{r_{c},t})$ into Equations (15) and (16) gives

[TABLE]

These are the values we use for our projections.

A.2 Estimating the Error Correlations

The model errors are defined by

[TABLE]

We have experimented with eleven different ways of estimating the correlations of the errors, which are the elements $A_{rs}$ of the matrix $\boldsymbol{A}$. Let $e_{r,t}$ denote the model error of region $r$ at time $t$. Let $\boldsymbol{\tilde{A}}$ denote a matrix where each element $\tilde{A}_{rs}$ is the empirical correlation between $e_{r\cdot}$ and $e_{s\cdot}$ over all time periods $t$, namely

[TABLE]

Furthermore, $\boldsymbol{\tilde{A}}$ truncated at zero will be denoted by $\boldsymbol{\tilde{A}_{[\geq 0]}}$. If positive definiteness is assured, it is denoted by $\boldsymbol{\tilde{A}^{*}}$.

We considered the following methods for estimating the matrix $\boldsymbol{A}$. In the first seven methods, all the within-country correlations are taken to be equal. The estimator of $\boldsymbol{A}$ is denoted by $\hat{\boldsymbol{A}}$. In all cases, $\hat{A}_{r,r}=1$ for all $r$, so in what follows, $\hat{A}_{r,s}$ refers to the cases where $r\neq s$.

$\hat{A}_{r,s}=\text{mean}\{a_{ij}\in\boldsymbol{\tilde{A}}_{[\geq 0]}\text{ and }i\neq j\}$ for all $r\neq s$. 2. 2.

$\hat{A}_{r,s}=\text{median}\{a_{ij}\in\boldsymbol{\tilde{A}}_{[\geq 0]}\text{ and }i\neq j\}$ for all $r\neq s$. 3. 3.

$\hat{A}_{r,s}$ is the Bayesian posterior mean of the intraclass correlation coefficient. 4. 4.

$\hat{A}_{r,s}$ is the Bayesian posterior mode of the intraclass correlation coefficient. 5. 5.

Similar to 3., but with errors divided by $\sqrt{1/n\sum_{r,t}e_{r,t}^{2}}$ with $n$ being the number of available errors. 6. 6.

Similar to 4., but with errors divided by $\sqrt{1/n\sum_{r,t}e_{r,t}^{2}}$ with $n$ being the number of available errors. 7. 7.

$\hat{A}_{r,s}=B/C$ for all $r\neq s$, where

[TABLE]

$N_{b}$ and $N_{c}$ the number of terms in the corresponding sum that are not missing, and $R$ is the number of regions. 8. 8.

The estimator of $\boldsymbol{A}$ is an approximation to the elementwise posterior median with uniform prior $U[0,1]$, namely:

[TABLE]

where $\bar{T}$ is the average number of time periods per region. It is a weighted average of the prior mean and the data. Note that this is our chosen method.

We will now show that if $\boldsymbol{\tilde{A}^{*}}_{[\geq 0]}$ is positive definite, then $\boldsymbol{A}$ is also positive definite. We can write

[TABLE]

where $\boldsymbol{B}$ is positive definite and $\boldsymbol{J}$ is the matrix all of whose entries are 1.

Now $\hat{\boldsymbol{A}}$ is positive definite if and only if $x^{\prime}\hat{\boldsymbol{A}}x>0$ for all $x\neq 0$. Now

[TABLE]

The first term on the right-hand side of (24) is positive by definition, since $\boldsymbol{B}$ is positive definitive. The second term is non-negative:

[TABLE]

Thus (24) is positive and so $\hat{\boldsymbol{A}}$ is positive definite. 9. 9.

Similar to 8. but with elements of $\boldsymbol{\tilde{A}}$ computed as

[TABLE] 10. 10.

Bayesian method introduced by Fosdick and Raftery, (2014): First, standardize $e_{r,t}$ by dividing the errors by $\sqrt{1/n\sum_{r,t}e_{r,t}^{2}}$ with $n$ being the number of available errors. Then the elements of the estimated correlation matrix $\hat{\boldsymbol{A}}$ are given as

[TABLE]

where $SS_{r}=\sum_{T}e_{r,t}^{2}$, $SS_{s}=\sum_{T}e_{s,t}^{2}$, and $SS_{r,s}=\sum_{T}e_{r,t}e_{s,t}$. Note that we are summing only over those time periods for which both countries, $r$ and $s$, have errors available. 11. 11.

Similar to 10., but using $(T+1)$ instead of $T$ in both the nominator and the denominator. This corresponds to the arcsin prior in Fosdick and Raftery, (2014). Note that a version of this method was tested where the errors were not standardized, but it performed less well, producing smaller correlations.

Fosdick and Raftery, (2014) found that correlations between countries were quite different for high and low TFR values. In light of this, we estimated two separate correlation matrices, one for the cases where the country had overall TFR 5 or above, and the other when the TFR was below 5.

The estimated correlation matrices resulting from methods 1.-7. have the same value for all off-diagonal elements. The elements of matrices resulting from methods 8.-11. differ from one another. In the latter case, all non-defined elements are set to the mean of the off-diagonal elements.

A.3 Out of Sample Validation Measures

This section provides detailed definitions for our out of sample validation measures.

We denote by $C$ the number of countries in our dataset, by $R_{c}$ the number of regions for country $c$, by $R$ the total number of regions, so that $R=\sum_{c=1}^{C}R_{c}$, and by $T$ the number of time periods over which we validate. Furthermore, $f_{r_{c},t}$ denotes the observed TFR, and $\hat{f}_{r_{c},t}$ denotes the point projection of the TFR (median of the predictive distribution), respectively, for region $r$ of country $c$ at time $t$.

The mean absolute error, MAE, is given by

[TABLE]

The bias is given by

[TABLE]

The Continuous Ranked Probability Score (CRPS) is an overall measure of the quality of a probabilistic forecast. If we are predicting the quantity $X$ (here the TFR), and produce the predictive distribution $F$, and observe a value $x$, then the CRPS is defined by:

[TABLE]

where $X^{\prime}$ and $X^{\prime\prime}$ are independent copies of a random variable with the distribution $F$ (Gneiting and Raftery, (2007), Eq. 21). We average the resulting values of CRPS across observations. Note that for the persistence method, the first part of the equation is zero. It is not simple to calculate the expectation, $E_{F}$, under the distribution $F$ analytically, so we did it by simulation. To obtain the expectation $E_{F}$, we sampled 5,000 values at random from the distribution $F$ and took the average of the corresponding predictands.

So far we have compared the methods for the TFR for all regions and for the average TFR for a country, taking the unweighted average over all regions in the country (Tables 1 and 2). Here, we add a comparison of the various correlation methods discussed in Appendix A.2 for the average TFR (Tables 3 and 4). The first row shows results when no correlation is taken into account. The number in parentheses of the following rows corresponds to the numbering of the eleven methods in the Appendix. In our study, we used method 8.