# Speech-based Age and Gender Prediction with Transformers

Felix Burkhardt1, Johannes Wagner1, Hagen Wierstorf1, Florian Eyben1, Björn
Schuller,1,2,3

###### Abstract

We report on the curation of several publicly available datasets for age and
gender prediction. Furthermore, we present experiments to predict age and
gender with models based on a pre-trained wav2vec 2.0. Depending on the
dataset, we achieve an MAE between $7.1$ years and $10.8$ years for age, and
at least $91.1$% ACC for gender (female, male, child). Compared to a modelling
approach built on hand-crafted features, our proposed system shows an
improvement of $9$% UAR for age and $4$% UAR for gender. To make our findings
reproducible, we release the best performing model to the community as well as
the sample lists of the data splits.

##  1 Introduction

The automatic detection of speaker age and gender has many use cases in human
computer interaction, for example for dialogue adaption or market research. In
contrast to subjective phenomena such as emotional arousal, the age of a
person may be objectively determined, like for example body size, by an exact
measurement. But, just like emotional arousal, age is only one of many factors
that influence the acoustic speech signal [30], and typically not to be
predicted to the year.

We curated several publicly available datasets with respect to age labels and
used them to train several age models based on a wav2vec 2.0 architecture. We
experiment on in- and cross-domain prediction, multi-head vs single head
models and the number of transformer layers to be used. Finally we report the
performance, using the 2010 paralinguistic challenge winner as a baseline.

Age and gender prediction based on machine learning as such has been
investigated numerous times in the past decades. A problem in this regard is
the lack of benchmark datasets that could be used to compare approaches. There
are several publicly available age annotated datasets like the SpeechDat II
corpus111https://tinyurl.com/speechdatii, CommonVoice [28], aGender [20],
Timit [21], VoxCeleb2 [27], or the NIST test set [19], but the studies we
found used only some of these, or not comparable train-development-test
splits. In addition, the authors usually only report either regression or
classification results, and use different metrics such as mean average error
(MAE), accuracy (ACC), precision, recall, or unweighted average recall (UAR).

Figure 1: Distribution of speaker age (#samples) in the datasets for the three
splits (CommonVoice age in mid-decades).

During the 2010 Interspeech Paralinguistic Challenge [22], age classification
was one of the topics. [9] report on the aGender dataset by fusing the results
of ensemble classifiers trained on subgroups of a larger feature set and get
$42.4\%$ UAR on four age groups, the baseline being $46.2\%$ UAR. [13] use a
similar configuration with respect to classifiers and feature sets to fuse
acoustic and metadata for child speech detection. Early studies are also based
on the aGender dataset [4, 6], both using gaussian mixture models
(GMM)/support vector machines (SVM) meta classifiers. In [3], the best system
reaches an F1 value of $.54$ with a Linear Discriminative Analysis (LDA) on
Hidden Markov Models (HMM) s modelling phonemes on the SpeechDat II corpus. A
human evaluation on a subset of the data reaches $.61$ F1 value. [18] describe
joint gender and age estimation using ivectors and Support Vector Regression,
reaching 4.7 years MAE on the NIST SRE 2010 telephony test set. [15] and [31]
both describe joint gender and age classification with a convolutional neural
network (CNN) based on the CommonVoice dataset and report a recall of $76\%$
on six mixed gender-age groups and $74\%$ UAR on twelve mixed gender-age
groups, respectively. Comparing linear and logistic regression on ivectors
with CNN s using spectrograms, [16] report $.98$ F1 for gender classification
and $9.44$ years MAE for age regression. They use part of VoxCeleb2 as a
dataset with age ground truth labels estimated mainly from a Wikipedia lookup,
and these labels are the basis for the VoxCeleb2 data used in this paper, see
Section 2.1 for details. [19] propose a long short term memory (LSTM)
recurrent network and report MAE of 6.58 years on the NIST test set. [17]
classified age using a wav2vec 2.0 model on the Timit dataset and report 5.54
years and 6.49 years MAE for male and female speakers, respectively. [33]
reached a similar performance on Timit utilising a QuartzNet architecture,
pre-trained on CommonVoice and VoxCeleb2, and fine-tuned on a joint age and
sex prediction. With respect to gender prediction, authors usually refer to
the biological sex. [29] report on the aGender dataset three-class problem and
reach 85% accuracy on the test set with a Random Forest classifier and MFCC
features. [32] use a pre-trained ResNet 50 and fine-tune it for gender on a
balanced sub-set CommonVoice. They report a recall of $.958$ on two gender
groups on the test set of CommonVoice and a similar recall for cross-corpus
performance.

With the paper at hand we see the following contributions:

  * •

We evaluate a novel system to estimate age and gender with fine tuned
transformer models

  * •

We present curated sample sets for train, development and test splits of
publicly available data sets and make them available to the research community

  * •

We compare a combined age and gender model to models specialized on a single
task

  * •

We present in-domain and cross-corpus results to examine the generalisability
of the proposed system

  * •

We investigate how many transformer layers are actually needed to properly
model the tasks

  * •

We release our best performing model to the public

##  2 Datasets

In order to test our approach on publicly available datasets, we considered
the ones mentioned in the introductory Section 1. Most of them have drawbacks:
The Timit and NIST datasets mainly contain young adults, Mozilla Common Voice
is labelled with self-reported decades as age, and aGender is only available
in $8$ kHz telephone quality, VoxCeleb2 might contain samples from the same
speakers, but recorded in different years. An overview on the age
distributions per split can be seen in Figure 1 and the numbers of samples and
speakers per dataset and split in Table 1. All datasets are available for non-
commercial research and the file lists can be accessed in the GitHub
repository that accompanies this
paper.222https://github.com/audeering/w2v2-age-gender-how-to.

###  2.1 VoxCeleb2

[16] report on a dataset which is based on a self-collected table for
VoxCeleb2 speakers. Because the authors (and the github repository) do not
provide exact sample lists but only the speaker splits, we limited the number
of samples per speaker to 20 samples. The original number of samples per
speaker in the VoxCeleb2 dataset has a mean value of 220. We re-used the test
set and split the train set randomly into 10% development speakers and 90% for
training. The age distribution for train and test splits is shown in Figure 1.
As said, we followed the splits used in [16] which have a rather large portion
of test speakers. Although the age peak is still on the young side, it
reflects the world age distribution better than, for example, Timit. As many
samples are quite long ($>20$ s), we used voice activity detection (VAD) to
segment the samples.

Dataset | train | devel | test
---|---|---|---
VoxCeleb2 | 30300 (1515) | 3120 (156) | 22100 (1105)
CommonVoice | 1729 (118) | 186 (13) | 1110 (79)
Timit | 1570 (157) | 170 (17) | 380 (38)
aGender | 29553 (324) | 2974 (35) | 20549 (239)
Table 1: Overview of the datasets: #samples and #speakers (in parenthesis).

###  2.2 CommonVoice

The dataset is the result of a public data collection by the Mozilla
foundation [28]. We used the de-validated collection as a basis for our
datasets. Because the number of samples per speaker varies strongly (Mean: 55,
STD: 296.6 ), we limited to 20 samples per speaker. Because the age
distribution is quite imbalanced, we furthermore tried to age-balance the test
and training splits by selecting at most 20 speakers per age decade and
gender, and then chose at most 7 speakers per age-gender group as a test set
and the others for training. As a development set, we randomly used 10% of the
training speakers.

###  2.3 Timit

The well known Timit datasets [21] contains $16$ kHz recordings of 630
speakers of eight major dialects of American English, each reading ten
phonetically rich sentences. Again, we tried to age-balance the test and
training splits by first selecting at most 40 speakers per age decade and
gender and then using 5 speakers per age-gender group as a test set,
disregarding many of the speakers in their twenties. As a development set, we
randomly used 10% of the training speakers.

###  2.4 aGender

The aGender dataset [20] has been collected via telephone especially with age
classification as a focus. It has been used in the Interspeech 2010
Pralinguistic challenge[22]. The winning paper was [23] and they reached
$53.86$ UAR and $81.57$ UAR for the age (4 classes) and gender (3 classes)
classification for the development set respectively with a meta classification
based on SVM and GMM. Because we wanted to be able to compare our results with
the challenge, we used the official development set as our test set and took
10% of the training speakers as a development set (the test set of the aGender
dataset is secret and not accessible). It is the only dataset that has been
collected solely for the purpose of biological age prediction and therefore
does contain a substantial amount of children and already a balanced age
structure.

##  3 Architecture

Figure 2: Proposed architecture built on wav2vec 2.0.

Our proposed architecture is depicted in Figure 2. It is built on wav2vec 2.0
[25] with two custom heads to predict age and gender, respectively. We do not
start training from scratch, but use the weights from a pre-trained model. In
our experiments we rely on _wav2vec2-large-robust_
333https://huggingface.co/facebook/wav2vec2-large-robust, a model pre-trained
on read speech from Libri-Light (60k hours) and CommonVoice (600 hours), but
also noisy telephone speech from Fisher (2k hours) and Switchboard (300 hours)
[26]. We could show that models fine-tuned on this variant are generally more
robust against noise compared to models that have seen only clean speech
during the pre-training [24].

As input to the heads we use the pooled hidden states (average pooling) of the
last transformer layer. Each head consist of a fully connected layer of size
1024, a dropout layer, and a final projection layer. In case of age, we
project to a single value predicting the _age_ in range 0 to 1, where 1
corresponds to a hundred years. In case of gender, we project to three values
expressing the confidence for being _child_ , _female_ , and _male_. During
evaluation we decide in favor of the class with the highest value.

For fine-tuning on the downstream task, we use the ADAM optimiser with a fixed
learning rate of $1\mathrm{e}{-4}$. Depending on the task we use two different
loss functions: concordance correlation coefficient (CCC) loss for age, which
we define as a regression problem; cross entropy (CE) loss for gender, which
we define as a multi-class classification problem. For backpropagation we use
the average of the two losses. We run for $5$ epochs with a batch size of $64$
and keep the checkpoint with best performance on the development set. As
proposed in [wagner2022dawn] we freeze the CNN layers but fine-tune the
transformer ones. When using the term fine-tuning, we will henceforth refer to
this partial fine-tuning. These models are trained using a single random seed,
for which the performance is reported.

In total, the model has $317.5$M parameters. On a three second long input it
performs $53.8$G MAC operations, which took $34.2\pm 5.2$ ms when measured on
a NVIDIA RTX A4000 (100 repetitions). In 4.3 we will discuss how these numbers
can be reduced. As learning framework we use PyTorch444https://pytorch.org/
and rely heavily on the transformer library by HuggingFace [34].

##  4 Experiments

We will now report results on predicting age and gender of a speaker from her
or his voice. We treat gender detection as a classification task and report
results in terms of accuracy (ACC)or unweighted average recall (UAR)for the
three classes child, female, and male. In case of age, which we model as a
regression problem, we report concordance correlation coefficient (CCC), and,
for the aGender age groups: ACC or UAR.

###  4.1 Single vs combined model

Figure 3: Performance of a model trained on a single task, i. e., either age
(MAE in years/100) or gender, and a combined model simultaneously trained on
both entities. We see that accuracy stays more or less the same.

We compare the performance of a _combined_ model trained on both tasks
simultaneously to that of models trained on a _single_ task, i. e., either age
or gender. For the latter, we use the same architecture described in Section 3
but remove the other head.

Results are summarised in Figure 3 and it is not difficult to recognise that
the single and combined model perform almost identical. We can conclude that,
although the combined model does not benefit from the information of the other
channel, it is well able to learn both tasks at once. Since this (almost)
halves the resources needed to run two separate models, a combined model
should be preferred in a multi-task setup. Throughout the remaining of the
paper, we will report results for the combined architecture without explicit
mention.

###  4.2 Cross-corpus evaluation

Figure 4: In the cross-corpus condition, the model has not seen data from the
dataset it is evaluated on. Results are compared with a model trained in an
in-domain fashion. The cross-domain model performs poorly in detecting females
(cf. confusion matrices) and predicts speakers younger than they are (cf.
distribution plots) (age: MAE in years/100).

In cross-corpus evaluation, a model is presented with data from an unknown
source. To simulate such a situation, we train a model on a single dataset.
For our experiment, we choose aGender as it is the only dataset that contains
a considerable amount of child speech. In Figure 4, we report the performance
of this model as _cross-corpus_ on the remaining datasets, namely CommonVoice,
Timit, and VoxCeleb2. For comparison, we also include _in-domain_ results by
the model trained on all datasets.

Removing in-domain data from the training generally leads to a performance
drop on all datasets. On CommonVoice, the effect is yet small: plus one year
for age and minus one percentage point for gender. Whereas on Timit and
VoxCeleb2 it is quite significant: MAE increases by four to five years for age
and ACC decreases by more than ten percentage points for gender. In the lower
part of Figure 4, we visualise gender and age predictions aggregated over the
three datasets. The confusion matrices reveal that the cross-domain model has
problems in predicting females, while the distribution plots show that age is
generally underestimated.

We conclude more and diverse data is needed to build a robust age and gender
model, especially when using an upsampled $8$ kHz dataset as training.

###  4.3 Varying the number of layers

Figure 5: Results with different number of transformer layers. For gender, a
single transformer layer seems sufficient, whereas for age, six layers provide
a good trade-off between speed and accuracy (age: MAE in years/100).

In the experiments reported so far, we have used all 24 transformer layers.
Dropping some of top layers reduces the footprint of a model. However, too few
layers may degrade the ability of the model to properly learn a task. To
investigate the effect of reducing transformer layers, we run experiments with
a varying number of layers. As we see in Figure 5, results generally improve
with more layers, though the effect is smaller for gender than for age. In
fact, a single transformer layer seems sufficient to adequately model gender.
For age, we can observe a considerable performance drop for less than six
layers, while with more than six layers, there is only a marginal increase.

We conclude that using six transformer layers provides a good trade-off
between accuracy and speed. This reduces the number of parameters by a factor
of $3.5$ to $90.8$M and inference time by a factor of to $3$ to $12.4\pm
3.3$ms.

###  4.4 Comparison classical modelling approach

| Age | Gender | Combined | Training
---|---|---|---|---
System | 4-class | 3-class | 7-class |
_Development_
baseline [23] | .56 / .55 | .82 / .87 | .54 / .54 | aGender
wav2vec 2.0 | .60 / .60 |  .84 / .92 | .58 / .60 | aGender
wav2vec 2.0 |  .61 / .61 |  .86 / .91 |  .59 / .61 | All
_Test_
baseline [23] | .52 / .51 | .83 / .86 | _N/A_ | aGender
wav2vec 2.0 |  .60 / .57 |  .86 / .90 |  .57 / .56 | aGender
wav2vec 2.0 |  .61 / .57 |  .87 / .88 |  .57 / .56 | All
Table 2: Comparison with baseline system based on handcrafted features.
Results are reported in terms of UAR/ ACC on the development set (upper part)
and test set (lower part). In case of age, we map the predictions of our model
to the four classes child, youth, adult, and senior. In the combined task, age
and gender are jointly represented by seven classes.

Finally, we compare performance to a classic modelling approach based on hand-
crafted features. As a baseline, we choose the winner system of the 2010
Interspeech Paralinguistic Challenge [22]. It implements a combination of
Gaussian Mixture Models (GMM) and Support Vector Machines (SVM), followed by
linear Gaussian backends and logistic regression-based fusion, which uses as
input a large feature set of acoustic, prosodic, and voice quality features.
For more information see [23].

In Table 2, we compare the performance of the baseline system with our best
performing model (24 layers) when trained either only on aGender or on all
datasets introduced in Section 2. In case of age, we map the continuous
predictions of our model to the four classes child, youth, adulate, and senior
as proposed by the challenge organisers. In addition, we include results for
the combined age/gender task with seven classes used in the challenge [22].
Depending on the task, UAR and ACC of the baseline is improved by 4–10
percentage points. Using all datasets during training provides an additional,
yet small boost.

We can conclude that deep learning improves the accuracy compared to a classic
modelling based on manual feature engineering.

###  4.5 Age/gender prediction of emotional data

When testing the published model on the Berlin emotional database [35], the
MAE for all samples age prediction is 8.35 and the UAR for binary gender
prediction 96.04. When only the neutral samples are used, the MAE drops to
5.94 and the UAR to 100. Clearly, the acted emotional expression jeopardizes
the quality of our model, as it has been trained on non-emotional data.

##  5 Conclusions and Outlook

We performed experiments on age and gender prediction based on four datasets
and a fine tuned transformer architecture.A model trained on all data sets,
together with the test, train, and develpment splits, has been made public and
can be used as a baseline for other authors. We will continue to investigate
age detection by using other model architectures and perhaps combining them
with expert features. Especially speech data from children is sparse and we
will look for such data or try to synthesise data with generative models.

##  6 Acknowledgements

This research has been partly funded by the European EASIER (Intelligent
Automatic Sign Language Translation) project (Grant Agreement number:
101016982) as well as the European MARVEL (Multimodal Extreme Scale Data
Analytics for Smart Cities Environments) project (Grant Agreement ID: 975337).

##  7 References

## References

  * [1] F. Eyben, M. Wöllmer and B. Schuller  “openSMILE –- the Munich versatile and fast open-source audio feature extractor”  In _Proceedings of the 18 th ACM international conference on Multimedia_, 2010, pp. 1459–1462
  * [2] B. Schuller, S. Steidl, A. Batliner, J.K. Hirschberg, A. Baird, A. Elkins, Y. Zhang1, E. Coutinho and K. Evanini  “The INTERSPEECH 2016 Computational Paralinguistics Challenge: Deception, Sincerity & Native Language”  In _Proceedings of the 17 th Annual Conference of the International Speech Communication Association, INTERSPEECH 2016_, 2016
  * [3] F. Metze, J. Ajmera, R. Englert, U. Bub, F. Burkhardt, J. Stegmann, C. Müller, R. Huber, B. Andrassy, J.G. Bauer and B. Littel  In _ICASSP, IEEE International Conference on Acoustics, Speech and Signal Processing - Proceedings_ 4, 2007
  * [4] T. Bocklet, A. Maier, J.G. Bauer, F. Burkhardt and E. Nöth  “Age and gender recognition for telephone applications based on GMM supervectors and support vector machines”  In _ICASSP, IEEE International Conference on Acoustics, Speech and Signal Processing – Proceedings_ , 2008
  * [5] F. Burkhardt, M. Eckert, W. Johannsen and J. Stegmann  “A database of age and gender annotated telephone speech”  In _Proceedings of the 7 th International Conference on Language Resources and Evaluation, LREC 2010_, 2010
  * [6] M. Feld, F. Burkhardt and C. Müller  “Automatic speaker age and gender recognition in the car for tailoring dialog and mobile services”  In _Proceedings of the 11th Annual Conference of the International Speech Communication Association, INTERSPEECH 2010_ , 2010
  * [7] M. Brückl  “Altersbedingte Veränderungen der Stimme und Sprechweise von Frauen” 7, Mündliche Kommunikation  Berlin: Logos Verlag, 2011
  * [8] Joel Shor, Aren Jansen, Ronnie Maor, Oran Lang, Omry Tuval, Félix Chaumont Quitry, Marco Tagliasacchi, Ira Shavitt, Dotan Emanuel and Yinnon Haviv  “Towards Learning a Universal Non-Semantic Representation of Speech”  In _Proc. Interspeech 2020_ , 2020, pp. 140–144
  * [9] Florian Lingenfelser, Johannes Wagner, Thurid Vogt, Jonghwa Kim and Elisabeth André  “Age and gender classification from speech using decision level fusion and ensemble based techniques”  In _Proceedings of the 11th Annual Conference of the International Speech Communication Association, INTERSPEECH 2010_ , 2010
  * [10] Nitesh V. Chawla, Kevin W. Bowyer, Lawrence O. Hall and W. Kegelmeyer  “SMOTE: Synthetic minority over-sampling technique”  In _Journal of Artificial Intelligence Research_ , 2002
  * [11] Florian Eyben, Klaus R. Scherer, Bjorn W. Schuller, Johan Sundberg, Elisabeth Andre, Carlos Busso, Laurence Y. Devillers, Julien Epps, Petri Laukka, Shrikanth S. Narayanan and Khiet P. Truong  “The Geneva Minimalistic Acoustic Parameter Set (GeMAPS) for Voice Research and Affective Computing”  In _IEEE Transactions on Affective Computing_ , 2016
  * [12] M Brückl and F. Heuer  “irrNA: Coefficients of Interrater Reliability – Generalized for Randomly Incomplete Datasets” R package version 0.1.4, 2018  URL: https://CRAN.R-project.org/package=irrNA
  * [13] Denys Katerenchuk  “Age group classification with speech and metadata multimodality fusion”  In _15th Conference of the European Chapter of the Association for Computational Linguistics, EACL 2017 - Proceedings of Conference_ , 2017
  * [14] Christopher Oates, Andreas Triantafyllopoulos, Ingmar Steiner and Björn W Schuller  “Robust Speech Emotion Recognition Under Different Encoding Conditions.”  In _INTERSPEECH_ , 2019, pp. 3935–3939
  * [15] Héctor A Sánchez-Hevia, Roberto Gil-Pita, · Manuel Utrilla-Manso, Manuel Rosa-Zurera and Manuel Utrilla-Manso  “Age group classification and gender recognition from speech with temporal convolutional neural networks”  In _Multimedia Tools and Applications_ , 2022
  * [16] Khaled Hechmi, Trung Ngo Trong, Ville Hautamäki and Tomi Kinnunen  “Voxceleb Enrichment for Age and Gender Recognition”  In _2021 IEEE Automatic Speech Recognition and Understanding Workshop (ASRU)_ , 2021, pp. 687–693
  * [17] Tarun Gupta, Duc Tuan Truong, Tran The Anh and Chng Eng Siong  “Estimation of speaker age and height from speech signal using bi-encoder transformer mixture model”  In _Proceedings of the Annual Conference of the International Speech Communication Association, INTERSPEECH_ 2022-September International Speech Communication Association, 2022, pp. 1978–1982
  * [18] Seyed Omid Sadjadi, Sriram Ganapathy and Jason W. Pelecanos  “Speaker age estimation on conversational telephone speech using senone posterior based i-vectors”  In _ICASSP, IEEE International Conference on Acoustics, Speech and Signal Processing - Proceedings_ 2016-May Institute of ElectricalElectronics Engineers Inc., 2016, pp. 5040–5044
  * [19] Ruben Zazo, Phani Sankar Nidadavolu, Nanxin Chen, Joaquin Gonzalez-Rodriguez and Najim Dehak  “Age Estimation in Short Speech Utterances Based on LSTM Recurrent Neural Networks”  In _IEEE Access_ 6 Institute of ElectricalElectronics Engineers Inc., 2018, pp. 22524–22530
  * [20] Felix Burkhardt, Martin Eckert, Wiebke Johannsen and Joachim Stegmann  “A Database of Age and Gender Annotated Telephone Speech”  In _Proceedings of the Seventh International Conference on Language Resources and Evaluation (LREC’10)_ Valletta, Malta: European Language Resources Association (ELRA), 2010
  * [21] J.. Garofolo, L.. Lamel, W.. Fisher, J.. Fiscus, D.. Pallett and N.. Dahlgren  “DARPA TIMIT Acoustic Phonetic Continuous Speech Corpus CDROM”  NIST, 1993
  * [22] Björn Schuller, Stefan Steidl, Anton Batliner, Felix Burkhardt, Laurence Devillers, Christian Müller and Shrikanth Narayanan  “The INTERSPEECH 2010 paralinguistic challenge”  In _Proceedings of the 11th Annual Conference of the International Speech Communication Association, INTERSPEECH 2010_ International Speech Communication Association, 2010, pp. 2794–2797
  * [23] Marcel Kockmann, Lukáš Burget and Jan Černocký  “Brno University of Technology System for Interspeech 2010 Paralinguistic Challenge”  In _Proceedings of the 11th Annual Conference of the International Speech Communication Association (INTERSPEECH 2010)_ 2010.9 Makuhari, Chiba, JP: International Speech Communication Association, 2010, pp. 2822–2825
  * [24] Johannes Wagner, Andreas Triantafyllopoulos, Hagen Wierstorf, Maximilian Schmitt, Felix Burkhardt, Florian Eyben and Björn W Schuller  “Dawn of the Transformer Era in Speech Emotion Recognition: Closing the Valence Gap”  In _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2023, pp. 1–13
  * [25] Alexei Baevski, Yuhao Zhou, Abdelrahman Mohamed and Michael Auli  “wav2vec 2.0: A framework for self-supervised learning of speech representations”  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2020, pp. 12449–12460
  * [26] Wei-Ning Hsu, Anuroop Sriram, Alexei Baevski, Tatiana Likhomanenko, Qiantong Xu, Vineel Pratap, Jacob Kahn, Ann Lee, Ronan Collobert, Gabriel Synnaeve and Michael Auli  “Robust wav2vec 2.0: Analyzing domain shift in self-supervised pre-training”  In _arXiv preprint arXiv:2104.01027_ , 2021
  * [27] Arsha Nagrani, Joon Son Chung, Weidi Xie and Andrew Zisserman  “Voxceleb: Large-scale speaker verification in the wild”  In _Computer Speech & Language_ 60 Academic Press, 2020, pp. 101027
  * [28] Rosana Ardila, Megan Branson, Kelly Davis, Michael Henretty, Michael Kohler, Josh Meyer, Reuben Morais, Lindsay Saunders, Francis M. Tyers and Gregor Weber  “Common voice: A massively-multilingual speech corpus”  In _LREC 2020 - 12th International Conference on Language Resources and Evaluation, Conference Proceedings_ European Language Resources Association (ELRA), 2020, pp. 4218–4222
  * [29] Sarah Ita Levitan, Taniya Mishra and Srinivas Bangalore  “Automatic identification of gender from speech”  In _Proceedings of the International Conference on Speech Prosody_ 2016-January International Speech Communication Association, 2016, pp. 84–88
  * [30] Susanne Schötz  “Acoustic analysis of adult speaker age”  In _Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics)_ 4343 LNAI Springer Verlag, 2007, pp. 88–107
  * [31] Anvarjon Tursunov, Joon Yeon Choeh and Soonil Kwon  “Age and gender recognition using a convolutional neural network with a specially designed multi-attention module through speech spectrograms”  In _Sensors_ 21.17 MDPI, 2021, pp. 5892
  * [32] Abeer Ali Alnuaim, Mohammed Zakariah, Chitra Shashidhar, Wesam Atef Hatamleh, Hussam Tarazi, Prashant Kumar Shukla and Rajnish Ratna  “Speaker gender recognition based on deep neural networks and ResNet50”  In _Wireless Communications and Mobile Computing_ 2022 Hindawi Limited, 2022, pp. 1–13
  * [33] Damian Kwasny and Daria Hemmerling  “Joint gender and age estimation based on speech signals using x-vectors and transfer learning”  In _arXiv preprint arXiv:2012.01551_ , 2020
  * [34] Thomas Wolf, Lysandre Debut, Victor Sanh, Julien Chaumond, Clement Delangue, Anthony Moi, Pierric Cistac, Tim Rault, Remi Louf, Morgan Funtowicz, Joe Davison, Sam Shleifer, Patrick Platen, Clara Ma, Yacine Jernite, Julien Plu, Canwen Xu, Teven Le Scao, Sylvain Gugger, Mariama Drame, Quentin Lhoest and Alexander Rush  “Transformers: State-of-the-Art Natural Language Processing”  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing: System Demonstrations_ Online: Association for Computational Linguistics, 2020, pp. 38–45
  * [35] Felix Burkhardt, Astrid Paeschke, M. Rolfes, Walter Sendlmeier and Benjamin Weiss  “A database of German emotional speech”  In _9th European Conference on Speech Communication and Technology_ 5, 2005, pp. 1517–1520  DOI: 10.21437/Interspeech.2005-446


# Forecasting Fertility using Parametric Mixture Models

Jason Hilton1
Erengul Dodd1
Jonathan J. Forster2
Peter W.F Smith1
Jakub Bijak1
1 Centre for Population Change, University of Southampton
2 Department of Statistics, University of Warwick
Correspondence<EMAIL_ADDRESS>

###### Abstract

This paper sets out a forecasting method that employs a mixture of parametric
functions to capture the pattern of fertility with respect to age. The overall
level of cohort fertility is decomposed over the range of fertile ages using a
mixture of parametric density functions. The level of fertility and the
parameters describing the shape of the fertility curve are projected foward
using time series methods. The model is estimated within a Bayesian framework,
allowing predictive distributions of future fertility rates to be produced
that naturally incorporate both time series and parametric uncertainty. A
number of choices are possible for the precise form of the functions used in
the two-component mixtures. The performance of several model variants is
tested on data from four countries; England and Wales, the USA, Sweden and
France. The former two countries exhibit multi-modality in their fertility
rate curves as a function of age, while the latter two are largely uni-modal.
The models are estimated using Hamiltonian Monte Carlo and the stan software
package on data covering the period up to 2006, with the period 2007-2016 held
back for assessment purposes. Forecast performance is found to be comparable
to other models identified as producing accurate fertility forecasts in the
literature.

##  1 Introduction

Human fertility is a complicated, dynamic process that is influenced by a wide
range of biological, social and economic factors as diverse as changes in
social norms around female education and employment; the availability and cost
of childcare; and wider trends in wages, disposable income and employment (see
Balbo, Billari, and Mills 2013 for a discussion). This complexity makes the
task of producing forecasts of future number of births difficult.

However, forecasts of fertility are vitally important in anticipating demand
for a variety of products and services. In the short term, demand for
maternity care and nursery school places and the aggregate cost of child
benefit payments all depend on the size of birth cohorts. Over longer
horizons, future fertility is a crucial determinant of overall population size
and age structure, and consequently will affect national accounts and fiscal
sustainability (Office for Budget Responsibility 2018).

As a result, methodologies that accurately capture the degree of uncertainty
surrounding forecasts may help manage the risks surrounding inevitable
deviations from point forecasts of fertility (Bijak et al. 2015). A predictive
forecast distribution can be combined with information about the costs of
over- or under-prediction and the degree of risk aversion of the decision
maker to allow a course of action to be decided upon, following the principles
of statistical decision theory (Bijak et al. 2015).

In several developed countries, fertility rates in recent decades can be
observed to take on a bi-modal shape when viewed as a function of age
(Chandola, Coleman, and Hiorns 1999). This suggests the possibility that the
underlying population is heterogeneous with respect to their fertility
behaviour (ibid), which may reflect, for example, the differential opportunity
cost of childbearing to women of different educational groups (Van Bavel 2010;
Billari and Philipov 2005). This paper describes an approach to the
forecasting of fertility that captures this feature of contemporary fertility
patterns. The method focuses first on forecasting the overall level of
childbearing, measured by the average number of births to mothers of
particular cohorts (that is, mothers born in the same year). Following
Chandola, Coleman, and Hiorns (1999), this summary measure of fertility is
subsequently decomposed over the fertile age range using a mixture of smooth
parametric curves. Extending their approach to a forecasting context,
predictions of future values are produced by forecasting forward both the
parameters of these curves and the summary measure of the level of fertility
using time series methods. Posterior predictive distributions coherently
incorporate both the underlying time series stochasticity and uncertainty
about model parameters.

The approach advocated in this paper has the advantage of parsimony, with
relatively few effective parameters needed for each cohort, and furthermore
constrains predictions of future fertility curves to the space of shapes
described by two-component mixtures. A number of authors have found such
mixtures to be useful in modelling fertility curves in populations in the
developed world (Chandola, Coleman, and Hiorns 1999; Peristera and Kostaki
2007; Bermúdez et al. 2012), but this approach adds dependence between
successive cohorts within a Bayesian hierarchical framework.

This paper makes several specific novel contributions to the literature on
fertility forecasting:

  * •

The existing work on parametric mixture modelling of fertility is extended to
a forecasting context, and a strategy for identifying the parameters of the
two components is set out.

  * •

The efficacy of employing different combinations of parametric forms within
the mixture model is examined.

  * •

Existing works that forecasting parametric fertility models without mixtures
(e.g. Congdon (1990); Congdon (1993); Knudsen, McNown, and Rogers (1993);
Mazzuco and Scarpa (2015)) employ a two stage approach whereby parametric
curves are first fitted to demographic rates, and time series models
subsequently fitted to the parameters of such models. The method set out in
this paper fits the parametric curves and time series models simultaneously
within a Bayesian framework, accounting for uncertainty in both elements of
the model in predictions.

  * •

Births are modeled using a negative binomial likelihood with smooth age-
specific over-dispersion, allowing for a better fit across the most fertile
ages.

Existing work on mixture forecasting in the statistical literature (e.g. Wong
and Li 2000; Li and Wong 2001; Wood, Rosen, and Kohn 2011) differs subtly to
the approach adopted in this paper. Firstly, the cited works apply to
univariate time series, whereas in this application the data are multivariate;
at each time point we have a set of correlated observations. Secondly, this
paper does not consider that observations themselves are drawn from a mixture
distribution, so that the likelihood for each point is a weighted sum of
density functions, but that particular parameters in a hierarchical
specification are weighted sums of parametric functions.

The rest of this paper proceeds by first describing the data used for the
forecasting application in the remainder of the paper (Section 2). Next, in
Section 3, existing approaches to the forecasting of fertility are introduced,
and areas in which further developments may be fruitful are identified.
Section 4 presents the model used throughout the paper, while Section 5
describes the model estimation process, and provides illustrative results of
the application of the model to data from England & Wales. Section 6 assesses
the performance of mixtures of different components using data from a
selection of developed country. The final section (Section 7) provides a
discussion of the overall contribution of the paper and directions for further
work.

##  2 Data and notation

This section describes the data used in this applications and sets out some of
the framework and notation used in the rest of the paper. Many of the details
will be familiar to demographers, but is included for completeness.

In most applications, the target variable for fertility forecasting is a
vector of age-specific fertility rates. Observed data on the number of births
$b$ is classified by age of the new-born’s mother at last birthday $x$ and the
year of observation $t$. The age-specific fertility rate for a given calendar
year can be defined in terms of expectation of the associated random variable
$B$.

$\displaystyle\begin{split}B_{xt}\sim\text{Poisson}(\lambda_{xt})\\\
E[B_{x,t}]=\lambda_{x,t}\\\ f_{xt}=\frac{E[B_{xt}]}{R_{xt}},\end{split}$ (1)

where $R_{x,t}$ is the cumulative population exposure to risk at age $x$ in
year $t$, measured in person-years. This latter quantity is often approximated
by the size of the female population aged $x$ at the midpoint of year $t$.

However, the cohort of the mother (that is, the year in which they were born)
is considered to be particularly important in the analysis of fertility,
because it is believed that decision to bear a child is influenced by past
life experiences and by the number of children a woman already has (Ryder
1964). This suggests that we might also be interested in an alternative
specification of the fertility rate $f_{xc}$, with $c$ indicating cohort.

When we wish to examine changes by cohort, data classified by year and age is
not ideal, as those recorded as aged $x$ at the midpoint of year $t$ may have
been born in year $t-x$ or in year $t-x-1$. For this work, data on births and
exposure were obtained from the Human Fertility Database (HFD) (Human
Fertility Database 2018). The HFD decomposes both birth and exposure age-
period data into approximate age-period-cohort cells (Jasilioniene et al.
2015), known as Lexis triangles because of the shape they make on the ‘Lexis
diagram’, a common tool for analysing demographic data (e.g. Minton et al.
2017). This decomposition allows data classified by cohort and age or by
cohort and period to be reconstituted from these lower-level building blocks
(Jasilioniene et al. 2015). Figure 1 is a Lexis diagram illustrating these
differing classifications. The life-course of an individual can be represented
on this diagram as a line which begins at the point along the horizontal axis
corresponding to their birth date, and is projected at an angle of 45 degrees
clockwise from the vertical. The points at which the line crosses vertical and
horizontal grid-lines represent transitions between calendar years and
birthdays respectively. The red square covers demographic events happening
during 2001 to the population aged 2 at their last birthday, and so
corresponds to Age-Period data. In contrast, the green parallelogram
corresponds to events happening to those born in during 2001 between the ages
of 1 and 2 (Age-Cohort data). Finally, the blue triangle is the ‘Lexis
Triangle’ corresponding to the same cohort for events occurring while they
were aged 2 during 2004 (Age-Period-Cohort data).

Figure 1: Lexis diagram illustrating cells corresponding to Age-Period, Age-
Cohort, and Age-Period-Cohort classified data

Age-Cohort data will be used for fitting models of fertility and producing
forecasts throughout the remainder of the paper. It is worth noting that at a
particular point in time, data for more recent cohorts (represented as
diagonal slices through the Lexis diagram) are incomplete, in that we only
have observations for younger ages. The forecasting problem is therefore to
both complete the fertility schedule of existing cohorts and to predict the
outcomes for new cohorts. It should be said that by using the HFD Age-Cohort
data as provided without accounting for the way it has been constructed from
Age-Period data, uncertainty in the final forecasts will be underestimated.
However, it is expected that this error will be small relative to other
sources of error.

Figure 2 shows fertility age profiles for selected cohorts for Sweden and for
England & Wales (with fertility aggregated for the latter two nations). In
both countries, the curve shifts to the right for later cohorts as women
postpone childbearing as the result of longer periods in education or greater
professional prospects (Billari et al. 2007). However, a notable difference
between the fertility curves for later cohorts for Sweden and England & Wales
is that in the latter case we see a distinct hump or shoulder for younger
ages, as identified by Chandola, Coleman, and Hiorns (1999). Such a feature is
also visible in Irish and to a lesser extent USA data (ibid). Fertility
projection models need to be flexible enough to capture this bulge at younger
ages in order to be able to adequately forecast future fertility rates for
these countries.

Figure 2: Cohort fertility schedules for England and Wales and for Sweden,
selected cohorts.

It is common to describe fertility for a given calendar year or cohort in
terms of summary measures. For cohort fertility, the most natural summary is
the area under the fertility hazard function for that cohort, which is the
average number of children born to a member of this birth cohort over the
course of their lifetime, assuming they survive to the highest possible age of
childbearing $X$ (Jasilioniene et al. 2015). This is generally termed
Completed Family Size (CFS), and can be approximated by the sum of the
individual age-specific rates:

$\displaystyle\begin{split}\theta_{c}&=\int_{12}^{50}\frac{\lambda_{ac}}{R_{ac}}da\\\
&\approx\sum_{x=12}^{49}f_{x,c}\end{split}$ (2)

The approximation is exact under the assumption the hazard and exposure are
constant within each age interval. The equivalent period-based measure is the
well-known Total Fertility Rate (TFR), which is constructed in the same way
from Age-Period data, substituting the $c$ cohort index for a $t$. This
measure approximates the average number of children born to a ‘synthetic’
cohort who experienced a given period’s age-specific fertility rates across
their whole fertile lifespan.

##  3 Modelling and Forecasting Fertility

A wide range of approaches have been deployed for modelling and forecasting
fertility, of which Booth (2006) and Bohk-Ewald, Li, and Myrskylä (2018a)
provide comprehensive reviews. An instructive area of the literature focuses
on merely modelling fertility rates rather than forecasting them. One approach
which has particularly inspired the current work was introduced by Hoem et al.
(1981), who used a range of density functions to approximate the shape of the
fertility curve for a given period. Density functions are useful in such a
context as they are parsimonious representations of uni-modal curves, and are
guaranteed to integrate to one. A minimum of two parameters are needed to
describe the location $\mu$ and scale $\tau$ of the curve, and a separate
parameter $\theta$ representing the TFR is employed multiplicatively to
determine the area under the curve and thus the overall level of fertility, so
that

$\displaystyle f_{x}=\theta\;g(x\;;\;\mu,\tau)$ (3)

where $g$ is a density function.

Chandola, Coleman, and Hiorns (1999) adapted this work to account for the
changes in observed data in the 1990s. In particular, the emergence of bi-
modal fertility curves in England and the United States necessitated the use
of two parametric components in order to obtain a reasonable approximation of
the data. Peristera and Kostaki (2007) and Bermúdez et al. (2012) went
further, including new families of density functions. Peristera and Kostaki
(2007) used split-normal distributions where the variance parameters take
different values either side of the mean, while Bermúdez et al. (2012)
employed Weibull density functions, which may be parameterised to have either
positive or negative skews.

Parametric approaches have also been employed for forecasting fertility, as
opposed to merely modelling its age structure, although these have tended to
involve only single-component models rather than mixtures. The double-
exponential nuptial model of Coale and McNeil (1972) was reparameterised and
adapted by Knudsen, McNown, and Rogers (1993) to allow for forecasting using
ARIMA models, while Congdon (1990) similarly employed the function of Hadwiger
(1940) in the forecasting of fertility for London boroughs. More recently, De
Iaco and Maggio (2016) used ARIMA methods to forecast forward parameters of a
gamma function fitted to Italian fertility, and furthermore used a Markov
field model to capture correlations in the error structure of this model. In
keeping with the approach adopted in this paper, Mazzuco and Scarpa (2015)
attempt to capture and forecast the bimodal structure of fertility using a
Flexible Generalisable Skew-Normal Distribution, although the parameters in
this model are hard to interpret demographically.

The debate in the demographic literature around whether fertility is better
understood as a cohort (Ryder 1964) or period (Ní Bhrolcháin 1992) phenomenon
has strongly influenced the forecasting literature. The argument revolves
around the observation that period fertility fluctuates in response to
economic or political circumstances in particular years, but often the
Completed Family Size (or Children Ever Born, CEB) of women in different
cohort remains fairly stable. This is because women adjust the timing of their
planned births (referred to as the tempo) in response to circumstances, but
the overall level (or quantum) is less subject to adjustment (Bongaarts and
Feeney 1998).

The greater stability of cohort fertility has led to a number of approaches
that focus on ‘completing’ cohort fertility. This involves inferring the
fertility rates for cohort which have already started childbearing. For
example, Evans (1986) completes cohort fertility based by predicting fertility
for ages 25-45 using regressions that use the level and approximate gradient
with respect to age of fertility rates between 15-25 as covariates. Similarly,
Congdon (1993) adapts his previous work (Congdon 1990) by using Hadwiger
functions to complete cohort fertility schedules from age 30 onwards, and
subsequently forecasts forward Hadwiger parameters using time series methods.
De Beer (1985) fits ARIMA models to smooth fertility rates in both the cohort
and age direction, so that information from both previous cohorts and from
earlier ages can be used to complete cohort fertility schedules. Cheng and Lin
(2010) use a more traditional Age-Period-Cohort model to incorporate effects
along both time-varying dimensions.

Myrskylä, Goldstein, and Cheng (2013) demonstrate a pragmatic method for
completing cohort fertility whereby for each age, the average change in the
fertility rate over the preceding 5-year period is assumed to persist for 5
years into the future. The authors describe this as the ‘freeze-slope’
approach, for obvious reasons, and this is contrasted with naive ‘freeze-rate’
approach, in which rates are held constant. Schmertmann et al. (2014)
incorporate elements of this approach within a conjugate normal-normal
Bayesian model that generates prior distributions that smooth fertility rates
in the both the age and cohort direction. The approach extracts three
principal components from a large collection of historical data taken from the
Human Fertility Database. One element of the constructed prior penalises
deviations from linear combinations of these principal components. An
additional element of the prior penalises deviations from both ‘freeze-rate’
and ‘freeze-slope’ projections. The priors are calibrated so that the size of
the penalty reflects the distribution of such deviations in historical data.
This approach is developed by Ellison, Forster, and Dodd (2018), who construct
a hierarchical Bayesian model maintaining the underlying assumptions of the
model of Schmertmann et al. (2014).

An alternative approach focuses more on extrapolation. Lee (1993) adapts the
principle component-based method of Lee and Carter (1992), originally
developed for mortality forecasting, to forecast fertility along the period
axis. Li and Wu (2003) extend this method using cohort as the primary axis,
while Hyndman and Ullah (2007) provide a more general approach within the
functional data framework. Wiśniowski et al. (2015) further adapt the Lee and
Carter (1992) approach as part of a comprehensive Bayesian population
forecasting framework.

The work of Ševčíková et al. (2015) is of particular interest because it has
been employed by the United Nations for the World Population Prospects since
2015. This method uses a Bayesian hierarchical model to describe the global
distribution of schedules of TFR evolution, noting that countries tend to
follow a common pattern of initially high fertility, decline, and later
stabilisation. Expert opinion is utilised to determine target patterns of
fertility over age (the so-called Proportionate Age Specific Rates) for each
country.

A recent comprehensive assessment of methods for completing cohort fertility
schedules established that few perform better than the naive ‘freeze rate’
technique (Bohk-Ewald, Li, and Myrskylä 2018a). Comparing the projected CFS
against observed data for each method over a wide range of countries and
forecast jump-off years, the authors found that only the methods developed by
Myrskylä, Goldstein, and Cheng (2013), Schmertmann et al. (2014), Ševčíková et
al. (2015), and De Beer (1985) consistently outperformed the naive forecast.
However, errors at the level of age-specific fertility rates were not
assessed, and performance in forecasting _new_ cohorts (which is essential for
population predictions) was also not tested.

The strategy adopted in this paper builds on the parametric approaches of Hoem
et al. (1981), Peristera and Kostaki (2007) and Bermúdez et al. (2012),
fitting explicit time series models for the parameters defining the model for
each year or cohort. Although the work of Bohk-Ewald, Li, and Myrskylä (2018a)
identified poor performance on the part of the parametric models of Hoem et
al. (1981) and Peristera and Kostaki (2007), the original models do not
involve any dependence between adjacent cohorts or calendar years. This work
rectifies this short–coming, while retaining the advantages of these models,
namely the parsimony of the parametric approach, the ability to separate the
overall level of fertility and the shape of the curve, and the restriction of
the shape of the fertility curve to a family of plausible values.

##  4 Model Description

This section describes in detail the different elements of the model, and the
justifications behind them. Starting with the definition of the likelihood,
births at each age are assumed to be distributed according to a negative
binomial distribution (parameterised according to the mean), so that

$\displaystyle B_{xc}\sim\text{Negative
Binomial}(R_{xc}f_{xc},\text{exp}(\phi_{x})),$ (4)

where $B_{xc}$ are the number of children born to women of cohort $c$ between
exact ages $x$ and $x+1$, $R_{xc}$ is the exposure to risk between these ages
for this cohort, $f_{xc}$ is the corresponding age-specific cohort fertility
rate, and $\phi_{x}$ is a parameter controlling the degree of over-dispersion
relative to the Poisson distribution, which varies with age $x$.

The fertility rate $f_{xc}$ is modelled as product of $\theta_{c}$, a
parameter describing the average number of births to cohort $c$ by maximum
childbearing age $X$, and $\xi_{xc}$, where the vector of
$\boldsymbol{\xi}_{c}$ parameters for each time period describes the age
distribution of the fertility rate schedule (termed the Proportionate Age-
Specific Fertility Rate (PASFR) by Ševčíková et al. (2015)), so that:

$\displaystyle\begin{split}f_{xc}=\;\theta_{c}\xi_{xc},\;\sum_{x=14}^{X}\xi_{xc}=\;1\;,\text{for
all}\;c\end{split}$ (5)

The age distribution $\boldsymbol{\xi_{c}}$ is modelled as a mixture of two
components, each parameterised by a location and a spread parameter. In
particular, as with Hoem et al. (1981), probability density functions are used
for the components as these have the convenient property that they integrate
to one. Alternatively, other functions with this property could also be used.
For a particular age $x$, an approximation of the mass of the distribution
between $x$ and $x+1$, weighted by the mixture parameter $\psi$, defines that
component’s contribution to $\xi$, so that

$\displaystyle\begin{split}\tau_{x,c}=&\psi_{c}\;g_{1}(x+\tfrac{1}{2};\;\mu^{(1)}_{c},\tau^{(1)}_{c})+(1-\psi_{c})\;g_{2}(x+\tfrac{1}{2};\;\mu^{(2)}_{c},\tau^{(2)}_{c})\\\
\psi\in&\;[0,1]\end{split}$ (6)

where $g_{i}(x;\mu,\tau)$ is the value at $x$ of the density function with
location $\mu_{c}^{(i)}$ and spread $\tau_{c}^{(i)}$, and $\psi_{c}$ is the
mixture parameter. The superscripts $(i)$ on the location and spread
parameters distinguish between the two mixture components.

###  4.1 Modelling evolution over time

The specification in Equation (6) leaves us with six unique parameters for
each time period: two location parameters, two scale parameters, a mixture-
weight parameter, and a parameter describing the overall level of fertility.
One practical problem is ensuring that the two components are distinct to
avoid identifiability and label switching issues. One strategy to circumvent
this is to force the location of one component be strictly greater than the
other. This can be achieved by re-parameterisation to use the sum and
difference of the locations rather than the raw values, combined with
enforcing constraints on these new parameters so that they are positive and
lie within reasonable bounds. Given the range of fertile ages (~14-50 year),
the sum of the locations should lie above $35$, and to allow identification of
the two mixture components, the gap is constrained to be at least $2$ years.

After this re-parameterisation, the resultant set of parameters
$\boldsymbol{\eta}_{c}$ can then by forecast forwards using auto-regressive
models,

$\displaystyle\begin{split}\boldsymbol{\eta}_{c}&=\\{\theta_{c},\;\psi_{c},\;\mu^{(s)}_{c},\mu^{(d)}_{c},\tau^{(1)}_{c},\tau^{(2)}_{c}\\}\\\
\boldsymbol{\eta}_{c}&\sim\text{N}\left(A\boldsymbol{\eta}_{c-1},\Sigma_{\eta}\right),\\\
\end{split}$ (7)

where superscripts $s$ and $d$ indicate the sum and difference of the location
parameters respectively. A special case of this class of models is the simple
independent random walk:

$\displaystyle\begin{split}\boldsymbol{\eta}_{c}=\boldsymbol{\eta}_{c-1}+\boldsymbol{\epsilon}_{c}\\\
\boldsymbol{\epsilon}_{c}\sim
N(\boldsymbol{0},\boldsymbol{\sigma^{T}_{\eta}}I).\end{split}$ (8)

Priors for the first elements of all the parameters modelled as time series
are chosen to be weakly informative. Flexibility as to the exact forecasting
model is possible; more complicated model classes, including ARIMA models
incorporating differencing and moving average elements, or stochastic
volatility models, could also be considered.

###  4.2 Model Components

Building on the models discussed previously by Hoem et al. (1981), Chandola,
Coleman, and Hiorns (1999), and Bermúdez et al. (2012), we examine three
possible parametric functions for the mixture components; the gamma density,
the Hadwiger or inverse Gaussian density, and the Weibull density. We wish to
adopt non-standard parameterisation for each of these functions so that they
are defined in terms of a location and a spread parameter. In this way,
elements of the model that are to be given time series priors have a
meaningful interpretation that is relatively consistent regardless of the
specific function used. For the gamma component, the mode and standard
deviation are used, while for the Hadwiger density the mean is instead adopted
as the location parameter. For the Weibull density, the form of the density
function necessitates other choices; the median and the distance from the
median to the upper quartile are instead used for the location and spread
parameters. The split normal distribution of Peristera and Kostaki (2007) is
not used, as although the authors find it fits the data relatively well, the
join in the function at the mode makes it conceptually unappealling.

Although Chandola, Coleman, and Hiorns (1999), Peristera and Kostaki (2007)
and Bermúdez et al. (2012) use the same density function for both the
components ($g_{1}$ and $g_{2}$) of their respective parametric mixtures,
there are advantages to combining different parametric functions within a
mixture, as the different functions used have different tail behaviour, and
may be better suited to approximating fertility at older or younger ages. For
example, the Weibull density function can skew in either direction, allowing a
steeper decline than is possible with the positively-skewed gamma and Hadwiger
distributions. These distinctions may help provide for better identification
of the two components. As a result, we examined all possible combinations of
the three functions.

###  4.3 Overdispersion

The Poisson distribution, a natural choice for modelling count data, has
variance equal to its mean. This assumption can be too restrictive in many
cases, and the use of the negative binomial distribution allows for over-
dispersion relative to the Poisson, which tends to lead to smoother estimates
of rates. However, the effect of over-dispersion is more pronounced in
absolute terms for larger counts, and so the use of a single over-dispersion
parameter for fertility can lead to a greater emphasis on the fit in the tails
of the age-specific rate distribution where counts are small than in the more
fertile stretches of the age range. To address this concern, we specify a
smooth, age-specific over-dispersion function using a penalised basis-spline
(Wood 2006), and allow the data to determine where greater over-dispersion is
needed. A random walk prior penalises deviations from smoothness in the basis
function coefficients (Lang and Brezger 2001):

$\displaystyle\begin{split}\phi_{x}&=\boldsymbol{\beta}\boldsymbol{B}(x)\\\
\beta_{i}&\sim\text{Normal}(\beta_{i-1},\;\sigma_{\phi}^{2}).\end{split}$ (9)

where $\boldsymbol{B}(x)$ is a matrix of B-spline basis functions.

Fertility data is often aggregated over open age intervals at the beginning
and end of the fertile age range. For example, in the UK, births to women of
ages 49 and above are reported in aggregate. In contrast to the Poisson case,
the sum of $n$ negative binomial distributions with different means does not
follow a standard distribution, so for simplicity the distribution of births
in these aggregate groups is approximated by a negative binomial distribution
with the appropriate mean and variance. Only a small proportion of births
occur to women in these age groups, so the approximation is unlikely to make a
significant difference to results.

##  5 Estimation and Results

The model has been tested on data from the Human Fertility Database for four
countries: USA; England & Wales (in combination); Sweden; and France. Each of
these countries have reasonably long time series available, although there are
distinct differences in the patterns evident in the data. The USA and England
& Wales display bimodality to a greater degree than France and Sweden.
Although the Human Fertility Database provides birth data from ages 12-55,
these data are extrapolated at the youngest and oldest ages where empirical
information is not available (Jasilioniene et al. 2015). As a result, HFD data
on births below age 15 and above 49 are aggregated for the purposes of model
fitting.

The latest 10 years of data were held back for each country and not used in
fitting, to provide for an out-of-sample assessment of the model described
above. The No-U-Turn Sampler (NUTS) variant of Hamiltonian Monte Carlo
(Hoffman and Gelman 2014), available in the Stan software package, was used to
obtained posterior samples of the model (Stan Development Team 2018).
Hamiltonian Monte Carlo simulates movement through the parameter space by
analogy to a physical system where the potential energy is equal to negative
log-posterior (Neal 2010).

Four chains, each consisting of 10000 iterations, were run for each
combination of model specification and country. Nine different models
(involving all combinations of the three parametric functions) were tested for
each of the four countries, giving a total of 36 different sets of posterior
samples. The first half of each chain was used as a warm-up period during
which the sampling parameters were allowed to adapt to the shape of the
posterior. With respect to the adaptation parameters, a slightly higher target
acceptance rate (0.95) and maximum NUTS tree-depth (12) was used than is
selected by default, providing for smaller integrator step-sizes and longer
NUTS trajectories (Stan Development Team 2019). The remaining samples, thinned
by a factor of 4 to prevent excessive memory use, were used for posterior
inference.

For the majority of model and country combinations, Gelman-Rubin split
$\hat{r}$ diagnostics are below the suggested 1.05 threshold and examination
of trace-plots indicate sampler convergence to the target distribution (Gelman
and Rubin 1992; Stan Development Team 2019). In a minority of cases (7 out of
36), diagnostics indicated problems were encountered during sampling. In
particular, sampling of models involving a Hadwiger distribution was found to
be difficult, particularly when the Hadwiger was used to model the component
corresponding to older ages.

In one case – the Weibull/Hadwiger model applied to French data – the two
parametric components were not distinguishable and obtaining posterior samples
was not possible. Multi-modal posteriors were obtained in 4 cases, with one or
more chains sampling from a mode located in a different area of the parameter
space to the others, generally because the mixture parameters approach a limit
for some cohorts with incomplete data. In these cases, either the mode
occupied by the majority of chains, or, where the four chains are divided
equally between two modes, the mode with the smallest Leave-One-Out
Information Criterion (LOOIC) (Vehtari, Gelman, and Gabry 2016) was chosen for
posterior inference, and the other chains discarded. Recomputed values of the
split $\hat{r}$ statistic indicate that in each case the remaining chains
appear to have converged to the selected mode. By choosing a mode based on
LOOIC for this small subset of models means that where posteriors are multi-
modal, only the mode that is best able to predict left-out data-points is
considered in the model comparison exercises that follow in Section 6.

Divergent transitions were encountered in two additional cases. Divergent
transitions indicate that the value of Hamiltonian in the No-U-Turn sampler
differs greatly over the course of a HMC trajectory due to integration errors
(Stan Development Team 2019). In practice this means that the curvature of the
posterior varies such that a single leap-frog integrator step-size is
insufficient to explore all areas of the posterior. In both cases, diagnostic
plots show that this problem was caused by significant posterior probability
being placed on low values of one of the time series innovation variance
parameters, which results in a narrower posterior in this area of the
parameter space than elsewhere. The sampler is unable to reach these areas of
the posterior, meaning the obtained posterior samples are truncated away from
zero for this parameter (Betancourt 2016). Re-parameterisation to a non-
centred model where standard normal innovations are sampled separately and
multiplied by the variance parameter may help resolve this problem (ibid), but
also makes enforcing constraints on the time series more difficult. Results
from these models are not reported, but qualitative examination of the
truncated posteriors suggest it is unlikely that these parameter combinations
would be among the best performing models.

The code used to fit the model and generate the results presented in this
section are available from github: https://github.com/jasonhilton/fert_mix

###  5.1 Illustrative Results for England & Wales

We first examine an example model fitted with data from the England & Wales.
In this case, we employ a gamma density for the parametric component
corresponding primarily to younger ages, and similarly a Weibull density for
older ages. The over-dispersion in births is allowed to vary by age in
accordance with Equation (9). All time-varying parameters are given random
walk priors. The posterior distributions are displayed in Figure 3. Birth data
is fully observed for all ages up until the 1966 cohort, who were 50 in the
last year of observation and thus deemed to have completed childbearing. The
posterior distributions for incomplete cohorts is informed by both their
partially complete childbearing experiences and the time series priors on the
model parameters, which penalise excessive deviations from past values. This
combination predicts a slight increase in completed family size in the short
term, before leveling off at around 2 children per woman.

Figure 3: Posterior Distribution of Model Parameters; gamma/Weibull Mixture;
England and Wales data

The mode of the lower gamma component fluctuates between ages 20 and 25, while
the median of the upper Weibull density is considerably higher, at around 33
years (Figure 3B). The weight given to the younger component decreases after
the parents of the baby-boom generation have completed their childbearing, as
individuals begin to delay fertility to later ages (Figure 3C). Uncertainty in
all parameters increases for later cohorts for which there are fewer or no
observations.

Figure 4 presents posterior distributions for the two weighted components and
their sum, the mean Age-Specific Fertility Rate, for six selected cohorts.
This plot shows how the size, position and weight of the two component changes
over cohorts. A slight lack of fit is observed for the 1950 cohort, which is
caused by the inability of the model to entirely capture the period
fluctuations associated with the baby bust. Qualitatively, it also shows the
model performs reasonably well at predicting the out-of-sample observations,
with most observations close to the posterior median, although a more formal
assessment taking into account the negative binomial uncertainty associated
with prediction of births is required to have confidence in its performance.
Similarly, the choices regarding the specification of functional components
should be justified in terms of predictive performance. The next section
examines the forecasting error and coverage of different variants of the model
to this end.

Figure 4: Posterior Distribution of Age-Specific Fertility Rates for Selected
Cohort. Vertical black lines for later cohorts represent forecast jump-off
points.

##  6 Assessing Model Variants

In order to choose between the different possible specifications of our model,
we wish to make a formal assessment of the performance of the model both in
predicting observations falling within the fitting period and in forecasting
new observations. To this end, for each possible configuration of functional
components (for example, gamma-Weibull or gamma-gamma) the model was estimated
on data up to 2006 for four countries separately.

The approximate Leave-One-Out (LOO) cross-validation approach introduced for
model assessment by Vehtari, Gelman, and Gabry (2016) and (implemented in the
loo R package) is used to assess the ability of the various model
specifications to predict data drawn from same data-generating process. This
method uses importance sampling to produce an approximation to the expected
log point-wise predictive density (ELPD) obtained by predicting every point in
the data-set using only the rest of the data, thus avoiding the need to refit
the model $N$ times leaving out each data point. Pareto smoothing is used
reduce the instability in the tails of the importance sampling weights
(Vehtari, Gelman, and Gabry (2016), termed Pareto Smoothed Importance
Sampling, or PSIS).

For some models, a small number of observations exhibited somewhat high values
of the k parameter of the Pareto distribution fitted to the importance
weights. This indicates the PSIS estimates of that points contribution to the
ELPD are unreliable, likely due to overly influential observations (Vehtari,
Gelman, and Gabry 2016). In these cases, the model was re-estimated a maximum
of three times with one or more unreliable observations left out, allowing the
contributions of these points to the ELPD value to be calculated directly and
combined with the existing PSIS estimates to provide the final ELPD value.
Where estimates for more than 3 data-points were unreliable, several points
were left out in each re-fitting to save on computation, with these points
grouped to maximise the distance between them within each refitting. In
practice, this process of refitting made little difference to the overall ELPD
scores or the model comparisons.

Figure 5 displays results in terms of the Leave-One-Out Information Criterion,
which transforms the ELPD so that it is on the same scale as the deviance and
other commonly used model assessment measures, such as the Deviance
Information Criterion and the Bayesian Information Criterion. The four sub-
plots each represent one of the four countries used in the assessment
exercise, while within each subplot, the first parametric component is
represented in the columns and the second components are in the rows. Thus,
the top-left cell refers to a model with a gamma component for younger ages
and a Weibull Component for older ages. Missing cells correspond to models for
which posterior samples are not available, as detailed above. The deeper blue
the cell of the plot, the better the model is at predicting left-out points
within the range of the data. For all countries, the gamma-Weibull model
obtains the best LOOIC value, although given the sampling variation of the
LOOIC metric, it is not always possible to definitively separate it from the
next best performing models. The Weibull-Weibull and gamma-gamma model also
perform relatively well across countries, although in the French case, the
gamma-Weibull model was the clear winner.

Figure 5: Pareto Smoothed Importance Sampling appoximation to the Leave-One-
Out Information Criterion for all combinations of Parametric Components for
England and Wales, France, USA, and Sweden

Although specifying models that are able to capture the properties of the
observed data well enough to predict left-out points is important, this paper
is primarily focused on forecasting. Forecasts of fertility rates are compared
to actual observed forecasts for the years 2007-2016 (As the data used is
indexed by cohort, the held-back sample is not strictly speaking for specific
years, but for an index t = c + a, which will consist of data from two
calendar years). The root mean squared error (RMSE) and the empirical coverage
over all age-specific rates is calculated for each model to assess both the
adequacy of the point forecasts and the calibration of the forecast
uncertainty. With respect to coverage, for each forecast we calculate the
proportion of corresponding observations of age-specific rates that fall
within the specified probabilistic predictive internals (90% and 50% in this
instance). Because we are observing only one set of correlated outcomes for
each country (age-specific rates for 2007-2016), rather than many independent
repetitions, this metric does not correspond exactly to the frequentist
concept of coverage, but it does however give an indication as to whether the
uncertainty in forecasts appear well calibrated compared to actual observed
deviations from predictions.

Table 1 shows a comparison for each model of the RMSE of the mean posterior
forecast over the 10 years of hold-back data 2007-2016. As well as variants of
the model set out in this paper, results are also reported for the naive
forecast, and for models identified by Bohk-Ewald, Li, and Myrskylä (2018b) as
performing well in predicting cohort fertility. The best performing models for
each country are in bold, while an NA indicates that obtaining results for
that combination of components was not possible. It is clear that parametric
models are competitive with the comparison models. The Weibull-Weibull model,
for instance, has lower forecast error than both the Myrskylä, Goldstein, and
Cheng (2013) and Schmertmann et al. (2014) models across all four countries,
although it is less successful than the freeze rate model for USA and Sweden.
The gamma-Weibull model, which was consistently amongst the models with the
strongest predictive power for points covered by the fitting period, similarly
outperformed all competitors in forecasting England & Wales and USA data.

Table 1: Comparison of the 10-year out-of sample forecast Root Mean Squared Error over different countries and models Model | RMSE
---|---
Model | E and W | France | Sweden | USA
Parametric Models
Gamma / Gamma | 0.0069 | 0.0049 | 0.0069 | 0.0090
Gamma / Hadwiger | 0.0083 | 0.0044 | 0.0068 | 0.0091
Gamma / Weibull | 0.0053 | 0.0071 | 0.0073 | 0.0092
Hadwiger / Gamma | 0.0142 | 0.0158 | 0.0233 | 0.0167
Hadwiger / Hadwiger | 0.0119 | NA | 0.0180 | 0.0160
Hadwiger / Weibull | 0.0059 | 0.0096 | 0.0089 | NA
Weibull / Gamma | 0.0075 | 0.0038 | 0.0048 | 0.0114
Weibull / Hadwiger | 0.0076 | NA | 0.0113 | 0.0110
Weibull / Weibull | 0.0052 | 0.0041 | 0.0067 | 0.0111
Other Models
Schmertmann et al. 2014 | 0.0073 | 0.0066 | 0.0082 | 0.0159
Myrskylä et al. 2013 | 0.0067 | 0.0063 | 0.0101 | 0.0121
Freeze Rates | 0.0068 | 0.0051 | 0.0040 | 0.0119

Tables 2 and 3 set out the the proportion of held-out observations that fall
within the 50% and 90% predictive intervals respectively. As was also clear
from the point forecast errors, US data appears relatively difficult to
forecast, with considerable under-coverage at the 50% interval for all models.
In contrast, for UK data, coverage for most models falls reasonably close to
the ideal values. Coverage properties of models with Weibull or gamma as the
first component appear reasonably good, in contrast to those models which use
Hadwiger density functions for younger ages, which tend to under-cover. The
Schmertmann et al. (2014) model also has strong coverage properties, perhaps
because this model is calibrated based on historical data across multiple
countries, while the model of Myrskylä, Goldstein, and Cheng (2013) appears to
systematically undercover, in agreement with the results in Bohk-Ewald, Li,
and Myrskylä (2018b).

Table 2: Comparison of the 10-year out-of-sample 50% empirical forecast coverage over different countries and models  | Model | E and W | France | Sweden | USA
---|---|---|---|---|---
Parametric Models
| Gamma / Gamma | 0.57 | 0.40 | 0.46 | 0.22
| Gamma / Hadwiger | 0.52 | 0.44 | 0.48 | 0.22
| Gamma / Weibull | 0.67 | 0.30 | 0.54 | 0.24
| Hadwiger / Gamma | 0.20 | 0.24 | 0.24 | 0.19
| Hadwiger / Hadwiger | 0.32 | NA | 0.32 | 0.19
| Hadwiger / Weibull | 0.50 | 0.48 | 0.62 | NA
| Weibull / Gamma | 0.56 | 0.61 | 0.63 | 0.37
| Weibull / Hadwiger | 0.52 | NA | 0.51 | 0.43
| Weibull / Weibull | 0.57 | 0.64 | 0.50 | 0.34
Other Models
| Schmertmann et al. 2014 | 0.49 | 0.60 | 0.60 | 0.36
| Myrskylä et al. 2013 | 0.18 | 0.08 | 0.11 | 0.11
Table 3: Comparison of the 10-year out-of sample forecast empirical 90% coverage over different countries and models  | Model | E and W | France | Sweden | USA
---|---|---|---|---|---
Parametric Models
| Gamma / Gamma | 0.80 | 0.81 | 0.77 | 0.77
| Gamma / Hadwiger | 0.80 | 0.90 | 0.77 | 0.79
| Gamma / Weibull | 0.95 | 0.71 | 0.92 | 0.49
| Hadwiger / Gamma | 0.58 | 0.62 | 0.61 | 0.52
| Hadwiger / Hadwiger | 0.82 | NA | 0.63 | 0.76
| Hadwiger / Weibull | 0.84 | 0.75 | 0.92 | NA
| Weibull / Gamma | 0.86 | 0.99 | 0.90 | 0.79
| Weibull / Hadwiger | 0.83 | NA | 0.64 | 0.79
| Weibull / Weibull | 0.82 | 0.99 | 0.82 | 0.84
Other Models
| Schmertmann et al. 2014 | 0.79 | 0.76 | 0.72 | 0.62
| Myrskylä et al. 2013 | 0.41 | 0.28 | 0.23 | 0.25

Interestingly, the model is able to capture fertility trends in Sweden and
France despite these countries displaying considerably more uni-modal
fertility schedules than in the USA or England & Wales. The ability of the
proposed family of models to capture differing shapes and trends in the tails
of the schedule may explain its utility in these cases.

Overall, for the countries studied it is not easy to state that one parametric
model specification is definitively better than another; the models best able
to predict left-out data from the period to 2006 did not necessarily have the
better RMSE or coverages over the period 2007-2016. However, the gamma-Weibull
model appears to offer a reasonable compromise, as it appears to consistently
produce good results across contexts.

##  7 Discussion

This paper has introduced a flexible family of models for the forecasting of
parametric mixtures, designed with the particular application of forecasting
bi-modal fertility in mind. The model can capture the shape of both bi- and
uni-modal age-specific fertility curve, and projects forward the level of
fertility together with parameters relating to the location, scale and mixture
weight of the two parametric components. The assessment of the proposed family
of models across a small range of countries show that they produce reasonable
forecasts of recent fertility trends that are competitive with models
identified in the literature as being best able to predict cohort fertility
(Bohk-Ewald, Li, and Myrskylä 2018b).

A relatively simplistic variant of the model has been proposed here, but there
is considerable potential for expansion, particularly in the area of the time
series models used for forecasting; fully-fledged ARIMA models will enable the
sort of immediate continuation of short-term trends that is a feature of the
models of Myrskylä, Goldstein, and Cheng (2013) and Schmertmann et al. (2014).
Capturing correlations between various model parameters may also enhance the
predictive capacities of the model, to the extent that such correlations
persist over time.

In order to more fully evaluate the wider applicability of the model to the
forecasting of fertility, a more comprehensive assessment exercise should be
carried out along the lines of Bohk-Ewald, Li, and Myrskylä (2018a), whereby
forecasts errors are evaluated over a large range of countries and across many
time periods. The present discussion has focused on the ability of the model
to predict cohort trends in fertility, but forecasting along the period
dimension is equally possible (e.g. Congdon 1993). Looking further afield, the
generality of the model described means that it may also be appropriate in
other contexts in which bi-modal curves must be forecast into the future.

## Acknowledgement

This work was supported by the ESRC Centre for Population Change - phase II
(grant ES/K007394/1) at the University of Southampton. The use of the IRIDIS
High Performance Computing Facility, and associated support services at the
University of Southampton, in the completion of this work is also
acknowledged. All the views presented in this paper are those of the authors
only.

## References

Balbo, Nicoletta, Francesco C Billari, and Melinda Mills. 2013. “Fertility in
Advanced Societies: A Review of Research: La fécondité dans les sociétés
avancées: un examen des recherches.” _European Journal of Population = Revue
Europeenne de Demographie_ 29 (1): 1–38. doi:10.1007/s10680-012-9277-y.

Bermúdez, S., R. Blanquero, J.A. Hernández, and J. Planelles. 2012. “A new
parametric model for fitting fertility curves A new parametric model for
fitting fertility curves.” _Population Studies_ 66 (3): 297–310.

Betancourt, Michael. 2016. “Diagnosing Suboptimal Cotangent Disintegrations in
Hamiltonian Monte Carlo.”

Bijak, Jakub, Isabel Alberts, Juha Alho, John Bryant, Thomas Buettner, Jane
Falkingham, Jonathan J. Forster, et al. 2015. “Letter to the Editor:
Probabilistic population forecasts for informed decision making.” _Journal of
Official Statistics_ 31 (4): 537–44. doi:10.1037/emo0000122.Do.

Billari, Francesco C., and Dimiter Philipov. 2005. “Women’s education and
entry into a first union. A simultaneous-hazard comparative analysis of
Central and Eastern Europe.” _Vienna Yearbook of Population Research_ 1
(2004): 91–110. doi:10.1553/populationyearbook2004s91.

Billari, Francesco C., Hans Peter Kohler, Gunnar Andersson, and Hans
Lundström. 2007. “Approaching the limit: Long-term trends in late and very
late fertility.” _Population and Development Review_ 33 (1): 149–70.
doi:10.1111/j.1728-4457.2007.00162.x.

Bohk-Ewald, Christina, Peng Li, and Mikko Myrskylä. 2018a. “Forecast accuracy
hardly improves with method complexity when completing cohort fertility.”
_Proceedings of the National Academy of Sciences_ 115 (3): 630–81.
doi:10.1073/pnas.1722364115.

———. 2018b. “Supplemental Material for ’Forecast accuracy hardly improves with
method complexity when completing cohort fertility’.” _Proceedings of the
National Academy of Sciences_ 115 (3): 630–81.
www.pnas.org/cgi/doi/10.1073/pnas.1722364115.

Bongaarts, John, and Griffith Feeney. 1998. “On the Quantum and Tempo of
Fertility.” _Population and Development Review_ 24 (2): 271–91.

Booth, Heather. 2006. “Demographic Forecasting: 1980 to 2005 in review.”
_Interational Journal of Forecasting_ 22: 547–81.
http://www.sciencedirect.com/science/article/pii/S016920700600046X.

Chandola, T, D A Coleman, and R W Hiorns. 1999. “Recent European fertility
patterns: fitting curves to ’distorted’ distributions.” _Population Studies_
53 (3): 317–29. doi:10.1080/00324720308089.

Cheng, Roger, and Eric S. Lin. 2010. “Completing incomplete cohort fertility
schedules.” _Demographic Research_ 23: 223–56. doi:10.4054/DemRes.2010.23.9.

Coale, A. J., and D. R. McNeil. 1972. “The distribution by age of the
frequency of first marriage in a female cohort.” _Journal of the American
Statistical Association_ 67 (340): 743–49. doi:10.1080/01621459.1972.10481287.

Congdon, Peter. 1990. “Graduation of Fertility Schedules: An Analysis of
Fertility Patterns in London in the 1980s and an Application to Fertility
Forecasts.” _Regional Studies_ 24 (4): 311–26.
doi:10.1080/00343409012331346014.

———. 1993. “Statistical Graduation in Local Demographic Analysis and
Projection.” _Journal of the Royal Statistical Society (Series A)_ 156 (2):
237–70.

De Beer, J. 1985. “A time series model for cohort data.” _J Am Stat Assoc_ 80
(391): 525–30.
http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?cmd=Retrieve&db=PubMed&dopt=Citation&list_uids=12267306.

De Iaco, Sandra, and Sabrina Maggio. 2016. “A dynamic model for age-specific
fertility rates in Italy.” _Spatial Statistics_ 17\. Elsevier Ltd: 105–20.
doi:10.1016/j.spasta.2016.05.002.

Ellison, Joanne, Jonathan J Forster, and Erengul Dodd. 2018. “Forecasting of
cohort fertility under a hierarchical Bayesian approach.”
https://eprints.soton.ac.uk/cgi/eprintbypureuuid?uuid=dead073b-5d4c-47af-b4bc-13c9dfdea9e6.

Evans, M. D. R. 1986. “American Fertility Patterns : A Comparison of White and
Nonwhite Cohorts Born 1903-56.” _Population and Development Review_ 12 (2):
269–93.

Gelman, Andrew, and Donald B Rubin. 1992. “Inference from Iterative Simulation
Using Multiple Sequences.” _Statistical Science_ 7 (4): 457–72.

Hadwiger, H. 1940. “Eine analytische reproduktionsfunktion für biologische
gesamtheiten.” _Scandinavian Actuarial Journal_ 3-4: 101–13.
doi:10.1080/03461238.1940.10404802.

Hoem, Jan M., Dan Madsen, Jorgen Lovgreen Nielsen, Else-Marie Ohlsen, Hans
Oluf Hansen, and Bo Rennermalm. 1981. “Experiments in Modelling Recent Danish
Fertility Curves.” _Demography_ 18 (2): 231–44.

Hoffman, M D, and Andrew Gelman. 2014. “The no-U-turn sampler: Adaptively
setting path lengths in Hamiltonian Monte Carlo.” _Journal of Machine Learning
Research_ 15 (2008): 1–31. http://arxiv.org/abs/1111.4246.

Human Fertility Database. 2018. “Human Fertility Database.” Max Planck
Institute for Demographic Research (Germany); Vienna Institute of Demography
(Austria). www.humanfertility.org.

Hyndman, Rob J., and Md. Shahid Ullah. 2007. “Robust forecasting of mortality
and fertility rates: A functional data approach.” _Computational Statistics
and Data Analysis_ 51 (10): 4942–56. doi:10.1016/j.csda.2006.07.028.

Jasilioniene, A., D. A. Jdanov, T Sobotka, E. M. Andreev, K. Zeman, and V. M.
Shkolnikov. 2015. “Methods Protocol for the Human Mortality Database.” Human
Fertility Database. https://www.humanfertility.org/Docs/methods.pdf.

Knudsen, Christin, Robert McNown, and Andrei Rogers. 1993. “Knudsen, Christin
McNown, Robert Rogers, Andrei.” _Social Science Research_ 22: 1–23.

Lang, Stefan, and Andreas Brezger. 2001. “Bayesian P-Splines.” Insitut Für
Statistik Sonderforschungsbereich 386. Munich: Ludwig-Maximilians-Universität-
München. https://epub.ub.uni-muenchen.de/1617/.

Lee, R D. 1993. “Modeling and forecasting the time series of US fertility: age
distribution, range, and ultimate level.” _International Journal of
Forecasting_ 9 (2): 187–202. doi:10.1016/0169-2070(93)90004-7.

Lee, Ronald D, and Lawrence R Carter. 1992. “Modeling and Forecasting U.S
Mortality.” _Journal of the American Statistical Association_ 87 (419).
American Statistical Association: pp. 659–71.
http://www.jstor.org/stable/2290201.

Li, Nan, and Zheng Wu. 2003. “Forecasting cohort incomplete fertility: A
method and an application.” _Population Studies_ 57 (3): 303–20.
doi:10.1080/0032472032000137826.

Li, W. K., and C. S. Wong. 2001. “On a logistic mixture autoregressive model.”
_Biometrika_ 88 (3): 833–46.

Mazzuco, Stefano, and Bruno Scarpa. 2015. “Fitting Age-Specific Fertility
Rates By a Skew-Symmetric Probability Density Function.” _Journal of the Royal
Statistical Society (Series A)_ 178 (1): 187–203.

Minton, Jon, Richard Shaw, Mark A. Green, Laura Vanderbloemen, Frank Popham,
and Gerry McCartney. 2017. “Visualising and quantifying ’excess deaths’ in
Scotland compared with the rest of the UK and the rest of Western Europe.”
_Journal of Epidemiology and Community Health_ 71 (5): 461–67.
doi:10.1136/jech-2016-207379.

Myrskylä, M, Joshua R Goldstein, and Yen-Hsin Alice Cheng. 2013. “New Cohort
Fertility Forecasts for the Developed World: Rises, Falls, and Reversals.”
_Population and Development Review_ 39 (1): 31–56.
http://onlinelibrary.wiley.com/doi/10.1111/j.1728-4457.2013.00572.x/abstract.

Neal, Radford. 2010. “MCMC using Hamiltonian Dynamics.” In _Handbook of Markov
Chain Monte Carlo_ , edited by Steve Brooks, Andrew Gelman, Galin Jones, and
Xiao-Li Meng Meng. Chapman; Hall / CRC Press.

Ní Bhrolcháin, Maire. 1992. “Period Paramount? A Critique of the Cohort
Approach to Fertility.” _Population and Development Review_ 18 (4): 599.
doi:10.2307/1973757.

Office for Budget Responsibility. 2018. “Fiscal Sustainability Report.” July.
Office for Budget Responsibility. doi:10.1142/S021812749900095X.

Peristera, Paraskevi, and Anastasia Kostaki. 2007. “Modeling fertility in
modern populations.” _Demographic Research_ 16\. doi:10.4054/DemRes.2007.16.6.

Ryder, N B. 1964. “The Process of Demographic Translation.” _Demography_ 1
(1): 74–82.

Schmertmann, Carl, Emilio Zagheni, Joshua R. Goldstein, and Mikko Myrskylä.
2014. “Bayesian Forecasting of Cohort Fertility.” _Journal of the American
Statistical Association_ 109 (506): 500–513. doi:10.1080/01621459.2014.881738.

Stan Development Team. 2018. “The Stan Core Library.” http://mc-stan.org.

———. 2019. “Stan Modeling Language Users Guide and Reference Manual.”
http://mc-stan.org/index.html.

Ševčíková, Hana, Nan Li, Vladimíra Kantorová, Patrick Gerland, and Adrian E.
Raftery. 2015. “Age-Specific Mortality and Fertility Rates for Probabilistic
Population Projections.” doi:10.1007/978-3-319-26603-9_15.

Van Bavel, Jan. 2010. “Choice of study discipline and the postponement of
motherhood in europe: The impact of expected earnings, gender composition, and
family attitudes.” _Demography_ 47 (2): 439–58.

Vehtari, Aki, Andrew Gelman, and Jonah Gabry. 2016. “Practical Bayesian model
evaluation using leave-one-out cross-validation and WAIC.” _Statistics and
Computing_ 27 (5). Springer US: 1–20. doi:10.1007/s11222-016-9696-4.

Wiśniowski, Arkadiusz, Peter W. F. Smith, Jakub Bijak, James Raymer, and
Jonathan J. Forster. 2015. “Bayesian Population Forecasting: Extending the
Lee-Carter Method.” _Demography_ 52: 1035–59. doi:10.1007/s13524-015-0389-y.

Wong, Chun Shan, and Wai Keung Li. 2000. “On a mixture autoregressive model.”
_Journal of the Royal Statistical Society (Series B)_ 62 (1): 95–115.

Wood, Sally, Ori Rosen, and Robert Kohn. 2011. “Bayesian Mixtures of
Autoregressive Models.” _Journal of Computational and Graphical Statistics_ 20
(1): 174–95. doi:10.1198/jcgs.2010.09174.

Wood, Simon N. 2006. _Generalised Additive Models: An Introduction with R_.
Boca Raton: Chapman; Hall / CRC Press.