Forecasting Australian fertility by age, region, and birthplace

National and sub-national fertility patterns in terms of age, geographical distribution, and birthplace are of interest to demographers, epidemiologists, health care personnel, government planners, and decision-makers (see, e.g., Rayer et al., 2009; Raftery et al., 2012; Kalasa et al., 2021). Future fertility is a crucial component of overall population size and age structure and, therefore, will affect national population accounts and fiscal sustainability (Eastwood and Lipton, 1999). According to the mother’s age group, age-specific fertility rates are ratios of the annual births per 1,000 of the female estimated resident population in that age group. As a widely used measure of the age pattern of fertility in a community, age-specific fertility rates are important demographic metrics that ought to be accurately forecasted in practice.

Parametric and nonparametric fertility modeling and forecasting methods have been proposed (see, e.g., Booth and Tickle, 2008; Shang, 2012, for comprehensive reviews). Recently, Bohk-Ewald et al. (2018) assessed the forecast accuracy of 20 commonly used fertility forecasting methods. In a separate study, Gleditsch et al. (2021) reviewed various fertility projection approaches applied to European data. In the parametric school of thought, Rogers (1986); Thompson et al. (1989); Bell (1992), and Keilman and Pham (2000) suggested fitting parametric distribution, such as Gamma distribution, to annual age-specific fertility rates. Then, the number of ages is reduced to a relatively smaller number of Gamma curve parameters; forecast the Gamma curve parameters using time-series techniques to forecast age-specific fertility rates. In the nonparametric school of thought, Bozik and Bell (1987) and Bell (1988, 1992, 1997) used the principal component analysis to obtain forecasts of age-specific fertility rates with only a small number of orthonormal principal components. From the viewpoint of mean squared error, these principal components and their associated scores capture the primary mode of information in the original data matrix.

Many fertility modeling applications in the literature only consider aggregated fertility at the national level. The fertility rates observed at the national level can be further disaggregated by characteristics of Australian women, namely age, birthplace, and region of residence. The Australian Government’s Center for Population publishes annual forecasts of fertility at the national, state, and territory levels following assumptions that recent trends of observed births will continue for all age groups in the coming years (see,e.g., Centre for Population, 2021b; McDonald, 2020). Recently, academics studied fertility rates in capital cities and regional areas in Australia (see, e.g., Wilson, 2017; Wilson et al., 2020) and fertility rates of immigrants (see, e.g., Baffour et al., 2021). To the best of our knowledge, however, a comprehensive study of age-specific fertility rates taking consideration of geographical region of residence and birthplace of Australian females is previously lacking. This paper aims to fill this gap by providing national and sub-national fertility forecasts following hierarchical structures disaggregated by geographical locations and birthplace.

The birthplace information presents a critical aspect linking fertility with migration. For instance: in 2020, the Australian total fertility rate was 1.65 births per woman, 1.68 for Australian-born women, and 1.55 for overseas-born women. Thus, immigrant fertility had little impact on the national fertility rate. The total fertility rate for the major immigrant groups did not differ significantly from the national fertility rate: 1.68 for United Kingdom-born, 2.07 for New Zealand-born, 1.58 for India-born, and 1.12 for China-born women (Australian Bureau of Statistics, 2021b). The lower rate for the China-born women was explained entirely by their lower fertility in the age range 15-24, at which ages many Chinese women included in the population are international students living in Australia temporarily and thus unlikely to give birth. With the fertility forecasts by birthplace, we can study the temporal change in fertility within a migrant group and compare fertility between different migrant groups.

The region of residence at birth registration reveals the geographical variability in the fertility experience of Australian women. As one of the most urbanized countries in the world, Australia has over 67% of its total population living in eight capital cities on 30 September 2021 (Australian Bureau of Statistics, 2022b). Significant geographic inequalities between different areas of the country in terms of income, education levels, employment opportunities, and population structure result in the substantial spatial variation of fertility in Australia (Evans and Gray, 2018). Modeling and forecasting age-specific fertility rates of various regions identify fertility differences between metropolitan areas and rural or remote regions of the country.

A difficulty in regional fertility forecasting is that the forecasts at the regional levels are constrained to be sum to the forecasts at the national level. This problem gives rise to forecast reconciliation. Forecast reconciliation occurs in national account balancing (Sefton and Weale, 2009), in tourist demand (Hyndman et al., 2011), in mortality forecasting (Shang and Hyndman, 2017; Li and Tang, 2019), in annuity price forecasting (Shang and Haberman, 2017), in solar energy forecasting (Yagli et al., 2020), in electricity load and demand forecasting (Jeon et al., 2019; Nystrup et al., 2020; Ben Taieb et al., 2020), and in tourism demand forecasting (Kourentzes and Athanasopoulos, 2019), to name only a few. This paper contributes a nonparametric approach to producing local age-specific fertility forecasts that consistently add to the national total.

A characteristic of age-specific fertility rates is the smooth shape over ages each year. We incorporate smoothness into the modeling and forecasting techniques to improve the accuracy of short-term forecasts. Implementing a smoothing technique can better capture the underlying trend of fertility changes, mitigating the influence of missing values and measurement noise in many sub-national series. With these aims in mind, we consider the functional time series method of Hyndman and Ullah (2007) to forecast national and sub-national age-specific fertility rates. As a generalization of Lee and Carter’s (1992) method, the functional time series method views age-specific fertility rates from a functional perspective, where age is treated as a continuum. We apply three existing forecast reconciliation methods, namely bottom-up, optimal combination, and trace minimization (MinT) methods, to reconcile fertility forecasts and improve forecast accuracy. The bottom-up method involves forecasting each of the disaggregated series and then using simple aggregation to obtain forecasts for the aggregated series (Kahn, 1998). The method works well where the bottom-level series have a strong signal-to-noise ratio. The optimal combination method obtains forecasts independently at all levels of the group structure. Then, linear regression is used with a least-squares estimator to combine and reconcile the forecasts optimally. The MinT method incorporates the information from a full covariance matrix of forecast errors in obtaining a set of coherent forecasts. By minimizing the mean squared error of the coherent forecasts across all levels of disaggregation, the method performs well on empirical data.

The remainder of this paper is structured as follows. In Section 2, we describe the data set — Australian national and sub-national fertility rates. In Section 3, we introduce a functional time series forecasting method for producing point and interval forecasts and then revisit the reconciliation methods in Section 4. We evaluate and compare point and interval forecast accuracies between the independent and grouped time series forecasting methods in Sections 5.1 and 5.2, respectively. We consider the averaged root mean squared forecast error for evaluating the point forecast accuracy. The reconciled forecast methods generally produce the most accurate overall point forecasts, particularly the optimal combination with trace minimization. For evaluating the interval forecast accuracy, we consider the mean interval score. We found that the base and reconciled forecasts provide similar forecast accuracy, although the latter one eases forecast interpretability. Conclusions are presented in Section 6, along with some reflections on how the methodology can be extended further.

Following Shang and Yang (2021), we use a hierarchical model coupled with grouped multivariate functional principal component decomposition technique to forecast age-specific fertility rates. The method models multiple sets of fertility functions correlated at any particular level of the chosen hierarchy. The common features in the data groups are considered in the estimation of the long-run covariance function, which is the sum of covariance function and autocovariance functions at various lags. The population-specific trend is captured by the principal components corresponding to the individual fertility series. Finally, forecasts of the empirical principal component scores are produced for each population before transforming back to the age-specific fertility rate forecasts. We present details of this functional time series forecasting method in the remaining of this section.

Using administrative division and birthplace as disaggregation factors, we can organize Australia’s sub-national age-specific fertility rates into two group structures shown in Figure 4. Both structures comprise three levels, with parent nodes in the top level. The middle level consists of at least one child node. For example, in Figure 4a, the NSW Coast ($R_{2}$) region as a parent node is comprised of areas including Hunter ($A_{2}$), Illawarra ($A_{3}$), and Mid-North Coast ($A_{4}$) as children. When multiple children exist for any particular parent node, female populations represented by these children nodes are highly correlated spatially, either within Australia (administrative division hierarchy in Figure 4a) or around the globe (birthplace hierarchy in Figure 4b). Hence, fertility series of children nodes belonging to the same parent are collectively modeled via the grouped multivariate functional time series method (Chiou et al., 2014; Shang and Yang, 2021; Shang and Kearney, 2022). We briefly describe this method to produce forecasts for Australian fertility rates.

Define $\{f^{(j)}(x)\}_{j=1,\cdots,\omega}$ as a set of smoothed sub-national fertility functions, with $\omega$ representing the number of children nodes in consideration, and $\omega=1$ corresponding to the special case of univariate functional time series. All $f^{(j)}(x)$ under consideration are square-integrable functions defined over the same interval of age $\mathcal{I}=[15,45]$. Let $\bm{f}(x)=\left[f^{(1)}(x),\dots,f^{(\omega)}(x)\right]^{\top}$ with $\omega\geq 2$ be a vector in the Hilbert space. Let $\mu^{(j)}(x):=\mathbb{E}[f^{(j)}(x)]$ denote the mean function for the $l^{\text{th}}$ population. The autocovariance matrix of $\bm{f}(x)$ at lag $\ell$ can be easily obtained as

$\displaystyle\gamma_{\ell}(x,z)=\text{Cov}\left[\bm{f}_{0}(x),\bm{f}_{\ell}(z)\right],$

with its $(l,j)^{\text{th}}$ ($1\leq l,j\leq\omega$) element given by

$\displaystyle\gamma^{lj}_{\ell}(x,z):=\mathbb{E}\Big{\{}[f_{0}^{(l)}(x)-\mu^{(l)}(x)][f_{\ell}^{(j)}(z)-\mu^{(j)}(z)]\Big{\}}=\text{Cov}\left[f_{0}^{(l)}(x),f_{\ell}^{(j)}(z)\right],\quad x,z\in\mathcal{I}.$

We can then obtain a long-run covariance function of $\bm{f}(x)$ as

$\displaystyle C_{lj}(x,z)=\sum_{\ell=-\infty}^{\infty}\gamma_{\ell}(x,z).$

For the $j^{\text{th}}$ population, denote a vector of long-run covariance functions as $\bm{C}_{j}(x,z)$ = $\Big{[}C_{j1}(x,z)$, $\cdots,C_{j\omega}(x,z)\Big{]}^{\top}$. We can then define an integral operator $\mathcal{A}:\mathbb{H}\rightarrow\mathbb{H}$ with the covariance kernel $\bm{C}(x,z)=\{C_{lj}(x,z),1\leq l,j\leq\omega\}$ such that

$\displaystyle(\mathcal{A}\bm{f})(x)=\int_{\mathcal{I}}\bm{C}(x,z)\bm{f}(z)dz=\begin{pmatrix}\langle\bm{C}_{1}(x,z)\,,\bm{f}(z)\rangle\\ \vdots\\ \langle\bm{C}_{\omega}(x,z)\,,\bm{f}(z)\rangle\end{pmatrix},$

where $\langle\bm{C}_{j}(x,z)\,,\bm{f}(z)\rangle=\sum_{l=1}^{\omega}\langle C_{jl}(x,z)\,,f^{(l)}(z)\rangle$ and $\bm{f}\in\mathbb{H}$. In practice, an empirical long-run covariance can be estimated as

$\displaystyle\widehat{C}_{\mathcal{h}}(x,z)=\sum_{\ell=-n}^{n}W\left(\frac{\ell}{\mathcal{h}}\right)\widehat{\gamma}_{\ell}(x,z),$

where $W$ is symmetric weight function, $\mathcal{h}$ is a bandwidth parameter; and the estimator of $\gamma_{\ell}(x,z)$ is defined in the form of

$$\widehat{\gamma}_{\ell}(x,z)=\begin{cases}\frac{1}{n}\sum_{t=1}^{n-\ell}\left[\bm{f}_{t}(x)-\widehat{\bm{\mu}}(x)\right]\left[\bm{f}_{t+\ell}(z)-\widehat{\bm{\mu}}(z)\right],&\ell\geq 0;\\ \frac{1}{n}\sum_{t=1-\ell}^{n}\left[\bm{f}_{t}(x)-\widehat{\bm{\mu}}(x)\right]\left[\bm{f}_{t+\ell}(z)-\widehat{\bm{\mu}}(z)\right],&\ell<0.\end{cases}$$

The bandwidth parameter $\mathcal{h}$ in the kernel function can greatly affect the estimator’s performance in finite samples. Therefore, we prefer to select $\mathcal{h}$ via a data-driven approach, such as the “plug-in” algorithm proposed in Rice and Shang (2017).

By Mercer’s lemma, there exists orthonormal sequences $\{\bm{\phi}_{k}=[\phi_{k}^{(1)},\cdots,\phi_{k}^{(\omega)}]^{\top}\}_{k=1,2,\cdots}$ of continuous functions in $\mathbb{H}$, and a non-increasing sequence $\lambda_{k}$ of positive numbers, such that

$\displaystyle\bm{C}(x,z)=\sum_{k=1}^{\infty}\lambda_{k}\bm{\phi}_{k}(x)\bm{\phi}_{k}(z),$

with the $(l,j)^{\text{th}}$ element $C_{lj}(x,z)$ of $\bm{C}(x,z)$

$\displaystyle C_{lj}(x,z)=\sum_{k=1}^{\infty}\lambda_{k}\phi_{k}^{(l)}(x)\phi_{k}^{(j)}.$

By the Karhunen-Loève theorem, the stochastic process for the $j^{\text{th}}$ population can be expressed as

$\displaystyle f^{(j)}(x)$

$\displaystyle=\mu^{(j)}(x)+\sum^{\infty}_{k=1}\beta^{(j)}_{k}\phi_{k}^{(j)}(x),$

(2)

where $\beta^{(j)}_{k}$ is the $k^{\text{th}}$ principal component score given by the projection of $[f^{(j)}(x)-\mu^{(j)}(x)]$ in the direction of eigenfunction $\phi_{k}^{(j)}$, that is, $\beta^{(j)}_{k}=\langle f^{(j)}(x)-\mu^{(j)}(x)\,,\phi_{k}^{(j)}(x)\rangle$.

Applying the decomposition of (2) to a time series of smoothed fertility functions $\{\bm{f}_{1}(x),\cdots,\bm{f}_{n}(x)\}$ gives

	$\displaystyle f_{t}^{(j)}(x)$	$\displaystyle=\mu^{(j)}(x)+\sum^{\infty}_{k=1}\beta^{(j)}_{t,k}\phi_{k}^{(j)}(x)$
		$\displaystyle=\mu^{(j)}(x)+\sum_{k=1}^{K}\beta^{(j)}_{t,k}\phi^{(j)}_{k}(x)+e^{(j)}_{t}(x),$

where $\Big{\{}\phi_{1}^{(j)}(x),\dots,\phi_{K}^{(j)}(x)\Big{\}}$ is a set of orthogonal basis functions commonly known as functional principal components for the $l^{\text{th}}$ population, with $\Big{\{}\bm{\beta}_{1}^{(j)},\cdots,\bm{\beta}^{(j)}_{K}\Big{\}}$ their related principal component scores and $\bm{\beta}_{k}^{(j)}=\left[\beta_{1,k}^{(j)},\dots,\beta_{n,k}^{(j)}\right]^{\top}$ for $k=1,\ldots,K$; $e_{t}^{(j)}(x)$ denotes the model truncation error function with mean zero and finite variance for the $j^{\text{th}}$ population. We select $K$ as the minimum of leading principal components reaching 95% of total variance explained (Shang and Hyndman, 2017), such that

$$K=\operatorname*{\arg\!\min}_{K:K\geq 1}\left\{\sum^{K}_{k=1}\widehat{\lambda}_{k}\Bigg{/}\sum^{\infty}_{k=1}\widehat{\lambda}_{k}\mathds{1}_{\left\{\widehat{\lambda}_{k}>0\right\}}\geq 0.95\right\},$$

where $\mathds{1}\{\cdot\}$ represents the binary indicator function. There are alternative methods for determining the number of retained components, such as those of Chiou (2012); Yao et al. (2005); Rice and Silverman (1991); Hall and Vial (2006), but they are beyond the main focus of this paper.

Collectively modeling $\omega$ women populations requires truncating at the first $K^{\text{th}}$ functional principal components of all time series as

$$\bm{f}_{t}(x)=\bm{\Phi}(x)\bm{\beta}_{t},$$

where $\bm{\beta}_{t}=\left\{\beta_{t,1}^{(1)},\dots,\beta_{t,K}^{(1)},\beta_{t,1}^{(2)},\dots,\beta_{t,K}^{(2)},\dots,\beta_{t,1}^{(\omega)},\dots,\beta_{t,K}^{(\omega)}\right\}^{\top}$ is an $(\omega\times K)\times 1$ vector of principal component scores and

$$\bm{\Phi}(x)=\left(\begin{array}[]{ccccccccc}\phi_{1}^{(1)}(x)&\cdots&\phi_{K}^{(1)}(x)&0&\cdots&0&0&\cdots&0\\ 0&\cdots&0&\phi_{1}^{(2)}(x)&\cdots&\phi_{K}^{(2)}(x)&0&\cdots&0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ 0&\cdots&0&0&\cdots&0&\phi_{1}^{(\omega)}(x)&\cdots&\phi_{K}^{(\omega)}(x)\end{array}\right)$$

a $\omega\times(\omega\times K)$ matrix of the associated basis functions. With the empirically estimated $\widehat{\bm{\Phi}}(x)$, we can then make $h$-step-ahead point forecasts as

$\displaystyle\widehat{\bm{f}}_{n+h|n}=\mathbb{E}[\bm{f}_{n+h}(x)|\bm{f}_{1}(x),\cdots,\bm{f}_{n}(x);\widehat{\bm{\Phi}}(x)]=\widehat{\bm{\mu}}(x)+\widehat{\bm{\Phi}}(x)\widehat{\bm{\beta}}_{n+h|n,k},$

where the empirical mean function is defined as $\widehat{\bm{\mu}}(x)=\frac{1}{n}\sum_{t=1}^{n}\bm{f}_{t}(x)$. We adopt a univariate time series forecasting method of Hyndman and Shang (2009) to obtain the forecast principal component score $\widehat{\bm{\beta}}_{n+h|n,k}$ (see, also Shang and Hyndman, 2017; Shang and Yang, 2021). The advantage of univariate forecasting method is its ability to handle nonstationary time series.

In addition to point forecasts, we also compute interval forecasts to assess forecast uncertainty. Specifically, the method of Aue et al. (2015) is considered to construct pointwise prediction intervals as follows.

1. 1)

   Using all observed data, compute the empirical functional principal components $\widehat{\bm{\Phi}}(x)$ with their associated estimated principal component scores $\{\widehat{\bm{\beta}}_{1},\cdots,\widehat{\bm{\beta}}_{n}\}$, where $\widehat{\bm{\beta}}_{t}=\big{\{}\widehat{\bm{\beta}}_{t,1}^{(1)},\dots,\widehat{\bm{\beta}}_{t,K}^{(1)},$ $\dots,\widehat{\bm{\beta}}_{t,1}^{(\omega)},\dots,\widehat{\bm{\beta}}_{t,K}^{(\omega)}\big{\}}^{\top}$. The in-sample forecasts are then constructed as

   $$\widehat{\bm{f}}_{\xi+h}(x)=\widehat{\bm{\mu}}(x)+\widehat{\bm{\Phi}}(x)\widehat{\bm{\beta}}_{\xi+h,k},\qquad\xi\in\{K,\cdots,n-h\},$$

   where $\widehat{\bm{\beta}}_{\xi+h}$ are $h$-step-ahead forecasts produced by univariate time series models based on $\widehat{\bm{\beta}}_{\xi}$.

2. 2)

   With the in-sample point forecasts, we calculate the in-sample point forecasting errors

   $$\widehat{\bm{\epsilon}}_{\zeta}(x)=\bm{f}_{\xi+h}(x)-\widehat{\bm{f}}_{\xi+h}(x),$$

   where $\zeta\in\{1,2,\cdots,M\}$ and $M=n-h-K+1$.

3. 3)

   Follow Shang (2018), we also use the nonparametric bootstrap approach to calculate pointwise prediction intervals. Determine a $\pi_{\alpha}$ such that $\alpha\times 100\%$ of in-sample forecasting errors satisfy

   $$\pi_{\alpha}\times\gamma_{\text{lb}}(x_{i})\leq\widehat{\bm{\epsilon}}_{\zeta}(x_{i})\leq\pi_{\alpha}\times\gamma_{\text{ub}}(x_{i}),\qquad i=1,\dots,p.$$

   Then, the $h$-step-ahead pointwise prediction intervals are given as

   $$\pi_{\alpha}\times\gamma_{\text{lb}}(x_{i})\leq\bm{f}_{n+h}(x_{i})-\widehat{\bm{f}}_{n+h|n}(x_{i})\leq\pi_{\alpha}\times\gamma_{\text{ub}}(x_{i}),$$

where $i$ symbolizes the discretized data points.

The approaches described above can be applied to both grouped structures shown in Figure 4 to compute point and interval forecasts at the most disaggregated levels. In this way, the base forecasts may not add to the national-level forecasts. In Section 4, we will discuss forecast reconciliation methods used to address this issue.

We consider administrative division and birthplace as disaggregation factors for modeling sub-national fertility rates in Australia. For a particular $j$^th ($j\in\left\{\text{COB}_{1},\ldots,\text{COB}_{19},Z_{1},\ldots,Z_{10},\text{World}\right\}$) birthplace, denote all the base forecasts at the Area level in year $t$ as $\bm{A}_{j,t}=\left\{f_{\text{A1}\ast\text{j},t},\ldots,f_{\text{A47}\ast\text{j},t}\right\}^{\top}$. Following the grouping structure shown in Figure 4a, we can express $\bm{f}_{j,t}$, a vector containing all series at all administrative division levels, as a product of matrices given by

$$\underbrace{\left[\begin{array}[]{c}f_{\text{AUS}\ast\text{j},t}\\ f_{{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\text{R1}\ast\text{j},t}}\\ f_{{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\text{R2}\ast\text{j},t}}\\ \vdots\\ f_{{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\text{R11}\ast\text{j},t}}\\ f_{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\text{A1}\ast\text{j},t}}\\ f_{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\text{A2}\ast\text{j},t}}\\ \vdots\\ f_{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\text{A47}\ast\text{j},t}}\\ \end{array}\right]}_{\bm{f}_{j,t}}=\underbrace{\left[\begin{array}[]{cccccccc}\frac{\textit{E}_{\text{A1}\ast\text{j},t}}{\textit{E}_{\text{AUS}\ast\text{j},t}}&\frac{\textit{E}_{\text{A2}\ast\text{j},t}}{\textit{E}_{\text{AUS}\ast\text{j},t}}&\frac{\textit{E}_{\text{A3}\ast\text{j},t}}{\textit{E}_{\text{AUS}\ast\text{j},t}}&\frac{\textit{E}_{\text{A4}\ast\text{j},t}}{\textit{E}_{\text{AUS}\ast\text{j},t}}&\cdots&\frac{\textit{E}_{\text{A47}\ast\text{j},t}}{\textit{E}_{\text{AUS}\ast\text{j},t}}\\ {\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\frac{\textit{E}_{\text{A1}\ast\text{j},t}}{\textit{E}_{\text{R1}\ast\text{j},t}}}&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}0}&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}0}&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}0}&\cdots&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}0}\\ {\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}0}&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\frac{\textit{E}_{\text{A2}\ast\text{j},t}}{\textit{E}_{\text{R2}\ast\text{j},t}}}&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\frac{\textit{E}_{\text{A3}\ast\text{j},t}}{\textit{E}_{\text{R2}\ast\text{j},t}}}&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\frac{\textit{E}_{\text{A4}\ast\text{j},t}}{\textit{E}_{\text{R2}\ast\text{j},t}}}&\cdots&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}0}\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ {\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}0}&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}0}&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}0}&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}0}&\cdots&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\frac{\textit{E}_{\text{A47}\ast\text{j},t}}{\textit{E}_{\text{R11}\ast\text{j},t}}}\\ {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&\cdots&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\ {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&\cdots&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&\cdots&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1}\\ \end{array}\right]}_{\bm{S}_{j,t}}\underbrace{\left[\begin{array}[]{c}f_{\text{A1}\ast\text{j},t}\\ f_{\text{A2}\ast\text{j},t}\\ \vdots\\ f_{\text{A47}\ast\text{j},t}\\ \end{array}\right]}_{\bm{A}_{j,t}}$$

or $\bm{f}_{j,t}=\bm{S}_{j,t}\bm{A}_{j,t}$, where $\bm{S}_{j,t}$ is a summing matrix whose elements include ratios of population at risk in various administrative divisions. A general notation for population at risk is given by $\textit{E}_{\text{g}\ast j,t}$, where $g$ stands for a particular administrative division listed in Table LABEL:tab:2. For example, $\textit{E}_{\text{A1}\ast\text{j},t}$ is the number of women in Sydney ($\text{Area}_{1}$) born in the $j$^th birthplace in year $t$, while $\textit{E}_{\text{AUS}\ast\text{j},t}$ represents the overall population of female residents in Australia born in the same place.

The grouping structure shown in Figure 4b splits Australian fertility rates at the national level by birthplaces. For a particular $g$^th ($g\in\left\{A_{1},\ldots,A_{47},R_{1},\ldots R_{11},\text{Australia}\right\}$) administrative division, denote all the base forecasts at the COB level in year $t$ as $\bm{\text{COB}}_{t}=\left\{f_{\text{g}\ast\text{COB1},t},\ldots,f_{\text{g}\ast\text{COB19},t}\right\}^{\top}$. We can express $\bm{f}_{g,t}$, a vector containing all series at all birthplaces, into a product of matrices as

$$\underbrace{\left[\begin{array}[]{c}f_{\text{g}\ast\text{World},t}\\ f_{{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\text{g}\ast\text{Z1},t}}\\ f_{{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\text{g}\ast\text{Z2},t}}\\ \vdots\\ f_{{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\text{g}\ast\text{Z{10}},t}}\\ f_{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\text{g}\ast\text{COB1},t}}\\ f_{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\text{g}\ast\text{COB2},t}}\\ \vdots\\ f_{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\text{g}\ast\text{COB19},t}}\\ \end{array}\right]}_{\bm{f}_{g,t}}=\underbrace{\left[\begin{array}[]{cccccccc}\frac{\textit{E}_{\text{g}\ast\text{COB1},t}}{\textit{E}_{\text{g}\ast\text{World},t}}&\frac{\textit{E}_{\text{g}\ast\text{COB2},t}}{\textit{E}_{\text{g}\ast\text{World},t}}&\frac{\textit{E}_{\text{g}\ast\text{COB3},t}}{\textit{E}_{\text{g}\ast\text{World},t}}&\frac{\textit{E}_{\text{g}\ast\text{COB4},t}}{\textit{E}_{\text{g}\ast\text{World},t}}&\cdots&\frac{\textit{E}_{\text{g}\ast\text{COB19},t}}{\textit{E}_{\text{g}\ast\text{World},t}}\\ {\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\frac{\textit{E}_{\text{g}\ast\text{COB1},t}}{\textit{E}_{\text{g}\ast\text{Z1},t}}}&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\frac{\textit{E}_{\text{g}\ast\text{COB2},t}}{\textit{E}_{\text{g}\ast\text{Z1},t}}}&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\frac{\textit{E}_{\text{g}\ast\text{COB3},t}}{\textit{E}_{\text{g}\ast\text{Z1},t}}}&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}0}&\cdots&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}0}\\ {\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}0}&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}0}&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}0}&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\frac{\textit{E}_{\text{g}\ast\text{COB4},t}}{\textit{E}_{\text{g}\ast\text{Z2},t}}}&\cdots&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}0}\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ {\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}0}&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}0}&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}0}&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}0}&\cdots&{\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\frac{\textit{E}_{\text{g}\ast\text{COB19},t}}{\textit{E}_{\text{g}\ast\text{Z10},t}}}\\ {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&\cdots&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\ {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&\cdots&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&\cdots&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1}\\ \end{array}\right]}_{\bm{S}_{g,t}}\underbrace{\left[\begin{array}[]{c}f_{\text{g}\ast\text{COB1},t}\\ f_{\text{g}\ast\text{COB2},t}\\ \vdots\\ f_{\text{g}\ast\text{COB19},t}\\ \end{array}\right]}_{\bm{\text{COB}}_{g,t}}$$

or $\bm{f}_{g,t}=\bm{S}_{g,t}\bm{\text{COB}}_{g,t}$. The summing matrix $\bm{S}_{g,t}$ is designed in the same fashion as $\bm{S}_{j,t}$, except for fixing a particular administrative admission $g$.

To ease notations, without causing confusion we may use general notations $\bm{B}_{t}$ for base level forecasts ($\bm{A}_{j,t}$ or $\bm{\text{COB}}_{g,t}$), $\bm{S}_{t}$ for either summing matrices ($\bm{S}_{j,t}$ or $\bm{S}_{g,t}$), and $\bm{f}_{t}$ for the reconciled forecasts at all levels of disaggregation ($\bm{f}_{j,t}$ or $\bm{f}_{g,t}$), respectively. In Sections 4.1 and 4.2, we shall present two forecasts reconciliation methods using these general notations.

The grouped functional time series forecasting method is applied to Australian sub-national fertility rates following both structures, illustrated in Figure 4. We then reconcile the obtained base forecasts via the bottom-up and optimal combination methods. To assess model and parameter stabilities over time, we consider an expanding window analysis of considered time series models (see Zivot and Wang, 2006, Chapter 9 for details), which is carried out in the following steps:

1. 1)

   Determine a particular grouped structure to disaggregate Australian total fertility rates. Specifically, birthplace is considered fixed when the group structure in Figure 4a is used, whereas administrative division is deemed to be fixed when the group structure in Figure 4b is used.

2. 2)

   Using fertility rates in 1981–2001 as a training set, we initially produce one- to 10-step ahead point and interval forecasts before reconciliation via the base, bottom-up, and optimal combination methods.

3. 3)

   Increase the training set with one more year of observations (i.e., fertility rates in 1981–2002), produce one-to nine-step-ahead point and interval forecasts, and compute the reconciled forecasts.

4. 4)

   Iterate the process with the size of the training set increased by one year each time until reaching the end of the data period in 2011.

The process produces 10 one-step-ahead forecasts, 9 two-step-ahead forecasts, and so on, up to one 10-step-ahead forecast. We report evaluation results of point forecast accuracy and interval forecast accuracy in Sections 5.1 and 5.2, respectively.