CC BY
41
A Bayesian Prediction for the Total Fertility Rate of Afghanistan Using the Auto-regressive Integrated Moving
Average (ARIMA) Model
Sayed Rahmi Khuda Haqbin and Athar Ali Khan
•
Department of Statistics and Operations Research Aligarh Muslim University, Aligarh-202002, India [email protected], [email protected]
Abstract
In this article, a simple methodology to predict the total fertility rate of Afghanistan via a Bayesian statistical analysis method has been applied. R- statistical analysis tool is used for data analysis. To forecast, the "bayesforecast" package is needed. It is a substitute package in Rfor the "forecast" package in the traditional (frequentest) statistical method. The Bayesian data analysis using the specific case of the general auto-regressive integrated moving average model (ARIMA) is processed as follows; As the first step, the stationarity of the given data-set is assessed, the time series has been made stationary by taking differences. After fitting several models, as the most appropriate fitted model, the ARIMA (1, 2,1) model has been fitted to the data. The accuracy of the fitted model is examined, and thereafter, the developed model is analyzed. The posterior computation is done, using the Markov Chain Monte Carlo (MCMC) simulation method. The method ultimately focuses on drawing relevant inferences including the 16 years prediction, and the results are; in general, found to be satisfactory.
Keywords: Fertility, Total Fertility Rate, Age -Specific Fertility Rate, Auto-regressive Integrated Moving Average Model, Stationarity, Hamiltonian Monte Carlo Algorithm, Prediction.
1. Introduction
Demography is the "study of the human population" [1]. The size of families and children's numbers had a substantially decreased trend over the 20th century. Especially after 1960, the family size and the number of children (fertility) have been showing a broadly favorable pattern for sustainable development. These changes in the population resulted from a declining fertility rate. The economic development level of a nation could be related to its fertility. The developed countries have been showing a lower level of fertility trend, along with a better education system, a revolutionary change in urbanization, greater wealth, and other factors. On the other hand, undeveloped countries have been showing a higher fertility trend, the increased level of fertility associated with a desire for family, labor, and caregivers in old age, lack of education, unawareness of contraceptives, strict adherence to traditional religious beliefs, and lower rates of female employees. Population growth and dependency ratios are directly affected by fertility changes. Over the past century, fertility has been one of the most determinants, of the population growth rate. Infant and child survival will be increased to provide greater access to education and health services, especially for women. By improving the participation of women in the labor force and empowering and reducing the number of children, women have over a lifetime. Low levels of fertility also contributed by improving maternal health, reducing the mortality rate of children, poverty alleviation, and economic development. Our objective; in this article includes
the future prediction of the fertility of a population. Clearly, for this purpose, we must use some summary index of fertility. The Total Fertility Rate is a composite measure of the fertility of a population, it is the most significant component for Fertility projections.
1.1. Total Fertility Rate
The total fertility rate calculates how many children a cohort woman will have. so it represents the average Family Size of the population. The total fertility rate is the sum of the Age-Specific Fertility Rate (ASFR) over all ages of the childbearing period and is expressed as "per woman" i.e., the sum divided by 1000. The average number of children a woman would have over her lifetime of reproductive potential, which is between the ages of 15 and 49, is the population's total fertility rate.
Forecasting the fertility rate has a long history of being the interest of demographers. The total number of live births in Australia, which is significant for population growth has been forecasted by using an auto-regressive moving average model [2]. Another most important demo-graphical forecasted characteristic could be the mortality forecast of the U.S. through the application of a time series method [3].
Fitted a bi-variate auto-regressive model which acted as a transform function model for forecasting the total fertility rate and mean ages of childbearing in the U.S annual data [4]. Age-specific mortality in the United States is predicted for the long term, with confidence intervals, using time series methods from 1990 to 2065. [5]. In a different publication, a parameterized model used the age profiles of fertility to quickly and simply describe age-specific rates. Future vital rate profiles are projected using time series models of the parameters, which reflect the temporal patterns of the age profiles [6].
A mixed two methods; statistical time series and mathematical demography and applied to predict stochastic models of fertility rates and the mortality rates of the U.S. dataset then by using the random-matrix product theory forecasted various measures of demography [7]. In a well-done study, a range of non-stationary dynamic factor models were used to jointly estimate the total fertility rate changes within the rather homogeneous clusters of Southern European countries in order to investigate the viability of the forecast. [8]. In this research, three methods for forecasting demographic processes are used. The population forecasting for the demographic data from 1980 to 2005 is particularly intriguing. Exogenous variable structural modelling based on theory and expectation [9]. A recent work applied ARIMA model for the total fertility rate of Pakistan is similar work of the paper from a classical perspective [10]. One another work lately has been conducted by using the two-time series models (ARIMA and ARIMAX) for live birth forecasting in Nigeria, which is a measure of population changes [11]. Numerous methods for predicting fertility have been mentioned; the ones that have been discussed thus far were mostly developed for wealthy nations that have already passed the fertility transition stage and currently exhibit low fertility patterns. Another important consideration is that the developed methods are used from traditional viewpoints. Pre-transition (high fertility), the fertility transition, and post-transition are the three phases of the changes in total fertility rate in a Bayesian projection model for the total fertility rate of country-specific forecasts of all countries (low fertility). The model is based on the most recent data from the United Nations Population Division [12].
In the same way, we aim to model the TFR of Afghanistan via the ARIMA model and achieve it is 10 years forecast by harnessing the attributes of the Bayesian paradigm. The procedures of this article are as follows:
• A real dataset of the total fertility rate (TFR) of Afghanistan has been provided.
• By taking two time differences the non-stationary time series came into the form of stationary series.
• The given data set has been analyzed from the Bayesian perspective using the "Bayesforecast" package, in R statistical language and software.
• "bayesforecast" package in R is a package provides many functions for data analysis in the Bayesian prospective.
• "bayesforecast" package is a Stan-endded pack.
• Stan is probabilistic programming language, which is made easy the data analysis, when there is complex probabilistic model.
• Stan work by taking sample of the equivalent to the posterior distribution fo the model, by different algorithms.
• NUTS(MCMC) algorithm, which is Stan sampling algorithm is applied in this article for the sampling purpose of the posterior distributions.
• The model fitting procedure is done by the help of Bayesian model fitting criteria "WAIC" and "LOOIC".
• The ACF and PACF plot of the original data and differenced data has been given, it could be a proper judge property for model fitting.
• After all the fitted and proper model which is decided for the given data is; ARIMA (1,2,1) model.
2. Methods
In this section, very briefly the method which is applied for this study are expressed. The Total Fertility Rate data-set in Afghanistan has been given. The assigned object is Bayesian prediction of Total Fertility rate in Afghanistan at least for 16 future years. For further analysis, we have been applied Bayesian method of forecasting. The Bayesian substitute R forecast package (in classical statistics), is bayesforecast R package.The stationarity of the data tested, in a Bayesian standpoint, about the model decided. ARIMA(1,2,1) model selected for further analysis. After all analysis series, we reach to a satisfactory result.
2.1. ARIMA(p,d,q) Model
A time series Yt is said to follows an auto-regressive integrated moving average (ARIMA) model if the dth difference Zt = VdYt is stationary and follows (ARMA) process. If Zt ~ ARMA(p, q) model, we say that Yt is an ARIMA (p, d, q) process. Where the three parameters; "p" in the Auto-regressive terms represents the current values depending on it's p — previous values (lag values), while "q" in the moving average terms represents the current deviation from the mean depending on q's previous deviations and "d" is the number of differences, by taking, it makes a non-stationary time series in the format of stationary. When stationary is not an issue, then we can define an auto-regressive moving average (ARMA) model as follows;
Zt = et + 01 Z1 + 02 Z2 +......+ QpZt — 01 e\ — 02 e2 —......— 0qet (1)
p q
Zt = et + £ 0'Zf—i — £ 0jet—j (2)
i=1 j=1
Where; et is a white noise process, and 0 parameter which is the coefficient of the authoritative terms Zt—i, and 0j the parameter, which is the coefficient of the moving average et—j terms. Using the backshift operator, we can write this more succinctly, and 01,..., 0p are the auto-regressive parameters to be estimated, 01,., 0q are the moving average parameters to be estimated and e\,., et a series of unknown random errors (residuals) that are assumed to follow a normal distribution. The (ARMA) model can be simplifies when Box-Jenkins back-shift operator applied, simplified (ARMA) model as follows: The auto-regressive AR(p) model;
Zt - <piZt-i - <2Zt-2 -......- $pZf-p = et
Where;
<(B) = 1 - <B -.........- <pBp and 0(B) = 1 - 01B - 02B2 -.....- 0qBq are polynomials in B
of degrees of p and q respectively. For the process to be stationary the root of <(B) = 0 must lie outside the unit circle, if the root of 0(B) = 0 is also outside the unit circle, it is called the invertible, there is a unique model corresponding to the likelihood in that case. we can write Equation 1 in the form of;
<(B)Zt = 0(B)et (3)
The Equation 3 is the Backshift form of the ARMA (p, q) and the experience suggests in data analysis models with AR and MA components often fit better data than the AR and MA pure. at
the same method from Equation 2 we can write;
p q p q
(1 - E <iBl)Zt = (1 - E 0jBj)et ^ ( E <iBl)Zt = ( E 0jBj)et i=1 i=1 i=1 i=1
<p (B)Zt = 0q (B)et (4)
The above (ARMA) model can be extended and written using differences;
BZt = Zt-1
The difference operator V = 1 - B is;
BZt = (1 - B)Zt = Zt - Zt-1. Power in the Backshift operator is equal to the possible difference that a time series can take;
B2 Zt = B(BZt) = B(Z-) = Zt-2 for the second difference we can write;
V2Zt = (1 - B)(1 - B)Zt = (1 - 2B + B2)Zt = Zt - 2Zt-1 + Zt-2 eventually the "d" difference can be write;
d
Yt - E Yt = (1 - B)dZt (5)
k=1
"d"is the order of differencing, in that case. The formal ARIMA (p, d, q) model in the Backshif form is produced by replacing "Zt" in the ARMA model with the differences described above
[15].
<p(B)(1 - B)dZt = 0q(B)et (6)
The Equation 6 is the Backshift format of the ARIMA(p,d,q) model, with "d" difference.
3. Data and Model Specification
There is a yearly data-set, the total fertility rate of Afghanistan, from the year of 1983-2022 reported by the; (World bank open data, free and open access to global development data), the data-set is taken from the; https://www.macrotrends.net/countries/AFG/afghan/fertility-rate. For a better accuracy the total fertility rate data-set in Afghanistan for the year of 1983-2022 are presented as bellow;
Table 1: Afghanistan - Historical Total Fertility Rate form 1982 to 2022.
1982-1992 7.450 7.454 7.458 7.461 7.465 7.469 7.472 7.474 7.477 7.479 7.482
1993-2003 7.516 7.551 7.585 7.620 7.654 7.560 7.465 7.371 7.276 7.182 7.041
2004-2014 6.900 6.760 6.619 6.478 6.272 6.066 5.859 5.653 5.447 5.269 5.090
2015-2022 4.912 4.733 4.555 4.414 4.273 4.133 3.992
The allocated projection of Afghanistan's total fertility rate for the next 16 years is the objective of this article. Graphical illustration is the first step in characterising a data-set for the analysis technique; thus, as a time series data-set has been provided, the most appropriate plot for the given data-set is time plot, as shown in Figure 1. From a plotted data, the below tips could be found;
• Patterns
• Unusual observations
• Change over time
• Relationship between variables
1990 2000 2010 2020
Time
Figure 1: Displays the time series plot of the Total Fertility Rate of Afghanistan.
As we mentioned above, Figure 1 shows a constant pattern, from 1982 to 1998, a very short increase, afterwards, then the total fertility rate of the country has been showing a decreasing trend till 2022. Over all Figure 1 is visualized a decreasing trend, therefore, for fitting ARIMA model, testing the stationarity is the versatile condition. In a stationary time series, the statistical properties (mean, variance...) are independent of time, it does not depend on the time at which the series is observed. Time series with seasonality and trend patterns are not stationary, the trend and seasonality affect the value of the time series at different point of time. A series with the same pattern in any point of time is a stationary time series, (e.g: white noise series). Stationary time series will have not predictable pattern in the long- term [14].
3.1. Differencing to eliminate a trend
Differencing is the way to remove trends (non-stationarity) in a time series data. A difference operator is defined as;
Vyt = yt - yt-1, (7)
as well as, more broadly, for order d
Vdyt = (1 - B)dyt,
(8)
Where B denotes the backshift operator (i.e., Bky — t = ytkfork > 1) A random walk, for example, is one of the most basic and extensively used time series models, although it is not stationary. We may create a random walk model like follows:
yt = yt—1 + wt, wt ~ N(0,q).
(9)
When the difference operator is applied to Equation 9, it produces a time series of Gaussian white noise errors wt:
V(yt = yt-1 + wt) yt - yt-1 = yt-1 - yt-1 + wt yt - yt-1 = wt [16].
3.1.1 Making use of the 'diff()' function
There are many methods for making a non-stationary (time series with trend, seasonality, cyclic,... ) data-set in the form of stationary. For the total fertility rate data-set, we applied two times differencing to making the time series stationary. Let's use R's diff() function to remove the trend from the total fertility rate time series. We will give some factors for identifying the stationarity of time series data; the correlation between a variable and itself at various time lags is known as the auto-correlation function (ACF). With the linear dependency of other yt-i, yt-2,..., yt-k eliminated, the partial autocorrelation function (PACF) assesses the linear correlation between a series yt and a lagged version of itself yt + k. In this part, for the very important point, we will
Figure 2: Displays the observed and differenced Total Fertility Rate of Afghanistan.
consider the Figure 2, and assesses it for further understanding of the process. In Figure 2 in the first row and first column is the time plot of the observed data, and the rest two other plots are the correlogram of the observed total fertility rate of Afghanistan, ACF and PACF respectively. The things we are going to learn from the correlogram are listed in bellow;
• The ACF at lag 0 is plotted as a
• reference point because it defaults to 1 (i.e., the correlation of a time series with itself).
• The blue horizontal lines represent the CIs for approximately 95%; and
• The auto-correlation is really high.
Figure 2 shows that although the partial autocorrelation at lag-1 is quite strong (equaling the ACF at lag — 1), the other values at lags > 1 are comparatively low, in contrast to what we observed for the ACF. Indicies for the time lag are again real-valued in the PACF plot, but there is no value for lag-0 since it is impossible to eliminate any intermediate autocorrelation between t and t — k when k = 0, and as a result, the PACF does not exist at lag-0. Now the second row of the Figure 2 shows a differenced and stationary of the total fertility rate of Afghanistan, for model identifying as per ACF and PACF which are given the second row of the Figure 2 second and third columns, the ACF show auto-correlation at lag-1 and PACF shows also correlation to lag-1 and decreasing correlation toward zero. Therefore, we will go to the other model fitting criteria of Bayesian approaches.
3.2. Model Specification
In a next step, after visualizing the time series data, checking stationarity, and removing trend of the data, we will approximate a model for the given data for further statistical analysis. The model used ought to call for the fewest parameters, whose values may be calculated from the observable series. Everything should be made as standardized but not simpler, according to a quotation from Albert Einstein in Parzen (1982, p. 68). For fitting a model to the data, few more steps are very important.
3.2.1 Bayesian Model Fitting Criteria
There are various methods for assessing and contrasting Bayesian models. The posterior predictive checks can be used to evaluate the model's fit to the data. It can be instructive to assess the prediction accuracy of each model under consideration, compare them, and decide what to do next if all of the models under consideration have discrepancies with the data [17].
As the total fertility rate data set is in use,we can do a WAIC and LOOIC criteria test for model selection.
Watanable-Akaike Infromation Criteria (WAIC): It has also been suggested to utilise the log pointwise posterior predictive density, with the posterior variance of the likelihood being used to derive the effective number of parameters.
n »
lppd = £ log p(yi\d)p(9\y)dd. (10)
i=1 J
For more details; [18].
elppdwaic = lppd — pwaic (11)
nn
WAIC = £ logE[p(yi\d,y)]) + 2 £ varpost{logpiy,\0)). (12)
i=1 i=1
For more details;[19].
Where the posterior mean of the likelihood of the ith observation is E[p(yi\theta,y)]) and £n=1 varpost(logp(yi\0)) the variance of each variable in the log predictive density summed over n data points. Prior knowledge is taken into account by the WAIC, and the posterior distribution is used in a non-normalway[16].
Leave-One-Out Cross Validation: When using the leave-one-out information criterion (LOO-IC), N-1 observations are used as the validation sample, the process is repeated N times, resulting
in a different observation being predicted each time, and the prediction results are added up to produce an estimate of expected log predictive density (elpd) that closely resembles the outcomes that would be obtained by gathering new data and performing the validation. We put the calculations into a R package named "loo".
Table 2: Displays the WAIC and LOOIC.
ARIMA Models LOOIC WAIC
ARIMA(2,2,2) -153.5 -155.5
ARIMA(0,2,2) -155.8 -156.0
ARIMA(2,2,1) -155.0 -156.8
ARIMA(2,2,0) -156.0 -156.3
ARIMA(1,2,0) -156.0 -155.5
ARIMA(0,2,1) -156.6 -156.9
ARIMA(1,2,1) -157.0 -156.9
ARIMA(1,2,2) -154.8 -156.1
The WAIC and LOOIC are given for several ARIMA models, as per Table 2 ; ARIMA (1,2,1) Model has been showing the proper and best model among all other models, with lesser WAIC and LOOIC.
4. Bayesian Formulation
A Bayesian model is a parametric model as classical(or frequentest) model, but in Bayesian an addition prior probability distribution for the parameters of the model have to be defined. And the parameters of the model treat as a random variable rather to unknown constants. For a
Bayesian model, the components are listed as bellow:
• the data, denoted by Z.
• the parameters, denoted by the Greek letters.
• the distributions of the model, given by a specification f (Z0 or F(Z0
• the prior distribution, specification of f (0) or F(0) or the distribution of the 0 .
f is the probability distribution (cdf), F is a symbol which denotes the cumulative distribution density(cdf), 0 is the prior.
4.1. Bayes' theorem
As a starting point for the Bayesian analysis, we would like to bring the Bayes' Theorem which is the important formula for Bayesian. Bayes' theorem is based on the conditional probability distribution:
f (elZ) = A™ (13)
f(x) is the unconditional (or prior) pdf of x. The proportionality format of the Equation 13 is;
f (e | Z) a f (Z I e)f (0) (14)
In the Equation 14 the f(0 | Z) is the posterior; the posterior is proportional to the likelihood times to the prior. The model is emphasized that the proportionality of the posterior is specifically related to the 0. As the demonstrator is independent of 0 we easily ignore it.
4.2. Choosing the Prior Distribution of the Parameters
For the four parameters the location, scale, auto-regressive, and moving average parameters,, we must describe a four-dimensional joint prior distribution in the general case. However, it is common practise to assume the parameters' prior independence, which means that before examining the data, we are unaware of whether the parameters are positively or negatively connected. We can just define four distinct priors for the four parameters with independence. We investigate weakly informative prior distributions.
TFR - ARIMA(1,2,1) mu0 = e0 - t(0,2.5,7) o0 - t(0,1,7) ar1 - N(0,0.5) max - N(0,0.5)
we can write the proportional prior density of the defined parameters;
e2
m(e0) a (1 + )-4
A
n2(fr1) a exp 0.5 n3(91) a exp 0.5 MP0) a (1 + f )-4
4.3. Likelihood Function for the Model
Z : Z1, Z2,..., Zt-j. are the observations, from Equation 2, the conditional density of Zt over
Zt-1, Zt—2,
..., Zt-p is given by;
f (Zt | Zt-1, Zt-2......Zt-p; e0, fru Qj, o2) a (1)exp^(Zt-e°-^=1 friZt-iQjet-j)2 (15)
From Equation 15 the likelihood function of the given density, can be approximated by its conditional forms;
T-d
L(Z | e0, fr, 9j, o2) a ^ f (Zt | Zt-1, Zt-2......Zt-p; e0, fr, Qj), (16)
t=p+i
For more details; [17].
For simplification we can write;
L(Z | e0, <pu Qj, o2) a (-1 )(J^A exp-i E^1 (Zt-e0-E<=1 friZt-i-j QfH)2 (17)
where; fr = fr1, fr2, •••, frp, Qj = Q\,Q2, •••,Qq• If one has a sample size of T, and would like to estimate a ARMA(p,q), there are many ways to estimate the parameters of the model: For the two times differenced Total Fertility Rate series in Afghanistan, the Likelihood can be defined from Equation 17 as;
f {Z | e0, fr1, Q1, o2) a (02exp(-£) eT=-2(Zt-e0Zt-1-Q1 et-1)2 (18)
where Zt obviously Vyt [18].
4.4. Posterior Distribution Function of the Model
As, it has been discussed in Section 4.1, Bayesian approach is the implication of the Bayes'theorem for further details Section 4.1. The prior and likelihood for the given data-set (total fertility rate) have been defined, now by applying the Bayes' theorem we can reach to the posterior distribution
of the Bayesian distribution of the data-set, once the posterior defined, it will be easy to apply any probabilistic programming language for fitting the posterior in a model and estimate the parameters. From Section 4.1, the posterior; posterior a likelihood x prior log(posterior) = log(prior) + log(likelihood)
g2 e2 g.2
p(e0,0i,o2 | Z) a (1 + -7)-4exp-0.5exp-015 (1 + )-4x
( _L gXp(— 212 ) £T=22 (Zt —C0 —4>1 Zt—1—91 et—1)2 (19)
From Equation 19, it is obvious which eventually we reach to a joint joint posterior distribution of all parameters, we can find the marginal density function of each parameter and then generate sample. As we obtained a joint posterior distribution of all parameters, we can find the marginal density function of each parameter and then generate sample. In this article, which is the total fertility rate data-set, in a Bayesian approach, we are going to analyse the data-set and estimate the parameters. As we know the most important step for a Bayesian data analysis is; to create a joint posterior function of all the parameters, now we have the joint posterior function of the parameters. There are many computational tools for the Bayesian data analysis, as we will discuss in Section 5, we only run the Code in R, in "bayesforecast" package' fucntion stansarima() fucntion, which is a Bayesian forecast'package in R, based of Stan, and using the MCMC simulation method, the sampling algorithm in this function is NUTS(MCMC).
4.5. Markov Chain Monte Carlo (MCMC)
The modern method for approximating complex forms of the posterior distribution is the Markov Chain Monte Carlo (MCMC) simulation method. The concept is akin to treating the posterior distribution as a population and then drawing samples from it repeatedly. When you draw a large enough sample (say 1,000), the sample distribution should be extremely close to the population distribution, as you learnt in basic statistics. The samples drawn in the above comparison are correlated, thus if the first sample is high, the second one is more likely to be high as well. This is necessary since there is no direct way to take samples from the posterior distribution, which typically has a very complex structure; instead, we have various procedures that can lead us to the posterior indirectly. Correlation between samples is usually not a big deal, except that we need to draw extra samples to compensate. In Sections 5.1 and 5.2 are given the values of the parameters and in Section 5.2 plots are showing the trance plot and the sample distribution of 4,000 samples using MCMC. The shape of the MCMC sample distribution approximates the shape of the posterior distribution [19].
5. Bayesian Inference Using NUTS(MCMC)
NUTS is an extension of Hamiltonian Monte Carlo (HMC) which is a type of MCMC. NUTS can even sample from the model with a greater number of parameters, it is a strength of NUTS sampler algorithm. Stan is a very good probabilistic programming language, the only disadvantage is that, Stan cannot sample for discrete parameters. NUTS uses a recursive algorithm to build a set of likely candidate points that spans a wide swath of the target distribution, stopping automatically when it starts to double back and retrace its steps. Hamiltonian Monte Carlo is a Markov Chain Monte Carlo (MCMC) algorithm that avoids the random walk behavior and sensitivity to correlated parameters that plague many MCMC methods by taking a series of steps informed by first order gradient information [20].
We will describe how to give the data to the model created by the R package "bayesforecast" in this part. For forecasting, the "bayesforecast" package in R is a Stan end-packed package. We must use R to give the data as a certain type to the "bayesforecast" package's function for ARMA
model fitting in order to accomplish this and obtain the outcome. In "bayesforecast" package, the model can be estimated using the guidelines below:
stan_sarima( ts,
order = с(1, 0, 0), seasonal = c(0, 0, 0), xreg = NULL, period = 0, chains = 4, iter = 2000, warmup = floor(iter/2), adapt.delta = 0.9, tree.depth = 10, stepwise = TRUE, prior_mu0 = NULL, prior_sigma0 = NULL, prior_ar = NULL, prior_ma = NULL, prior_sar = NULL, prior_sma = NULL, prior_breg = NULL, series.name = NULL,
Using the function stan_sarima()'s from "bayesforecast" package in R the default setting are; ts =, The ARIMA(p,d,q) order, chains is equaling 4, iterations per chain equaling 1000, and 1000 warmups. used as a warmup to eliminate dependencies in the initial settings. We fit the selected model to the total fertility rate; by the function above.
5.1. The Estimated Results
The given Table 3: shows the numerical estimated values for the parameters of the fitted model ARIMA (1,2,1) on the total fertility rate of Afghanistan.
Table 3: Displays the estimated posterior distribution of the parameters
Mean Se %5 % 95 Ness Rhat
Intercept -0.0103 0.0001 -0.0215 -0.0007 4282.407 1.0000
0.0270 0.0001 0.0221 0.0329 3940.089 1.0015
Ф1 0.9209 0.0008 0.8339 0.9872 4034.682 1.0000
01 0.0658 0.0024 -0.1909 0.3134 4109.169 1.0005
loglik 77.7385 0.0323 73.8071 80.3939 3808.628 1.0034
Table 3 provides a summary of the statistics for the marginalized posterior distribution as well as the MCMC convergence diagnostics. In line 3, for instance, the marginalized posterior distribution for parameter 01 is presented. Now, let's examine each column. The name of the parameter is displayed in the first column. The posterior mean, also known as the mean of the marginalized posterior distribution, is the value in the second column. The arithmetic mean of all the draws in this case, a total of 4000 is used to calculate this. For instance, the posterior mean of 0i is 0.0658 , and this finding shows that, on average, the baseline annual income grows by 0.0658with each extra year of work experience. The Monte Carlo standard errors of the Mean are shown in the third column as Se. By dividing the SD by the square root of N — eff, it is calculated. The quantiles of the posterior marginalized distributions are displayed in columns
four and five. These are calculated using the draw's quantiles. Stan determines the effective number of drawings as Neff in column 6 using the autocorrelation of the draws. We believe this value should be at least 100 in order to estimate distributions and perform other statistical calculations. Additionally, a small value for this parameter suggests that it has probably not yet converged; this information could be used to improve the model. An indication of MCMC convergence, R — hat, is computed for each parameter and is shown in column seven. Typically, it compares the sample variance of each chain using findings from many chains. According to Section 11.4 [21].
The MCMC is said to have converged if it meets the condition Rhat < 1.1 for all the parameters. We must confirm that the MCMC is converged before evaluating the results. Beginners frequently skip over this stage in favour of moving on to the analysis' next step, when they could use the posterior mean, plot the histogram of draws, or quickly interpret the findings using an MCMC that hasn't even begun to converge. Without a doubt, avoid this. The try — and — error procedure must be repeated numerous times before the MCMC converges. Eventually the ARIMA (1,2,1) model for the time series data as per Equation 1 can be written:
Yt = -0.0103 + 0.92yt-1 + 0.66et-1 (20)
From Table 3, the simulation estimate for E(fr1 |y) is 0.92009, according the 90% credible interval which is (0.8339 , 0.9872) of the Table 3 statistically significant. in the same way, the simulation estimate for E(Q1 |y) is 0.06258, it is not statistically significant 90% probability interval for the variable(parameter) Q1 is (-0.1909 , 0.3134) . Intercept and the posterior standard deviation values are given respectively. In Bayesian statistics, the credible interval is analogous to the confidence interval, but with radically distinct implications. A 90% credible interval, for example, is an interval that has a 90% chance of containing the true value of the parameter, an interpretation that is frequently and incorrectly connected with the confidence interval. Because the population parameter is fixed, frequentists cannot use probability for the parameter; instead, only the sample is probabilistic, and a 90% confidence interval must be interpreted in the sense that 90% of the intervals constructed with repeated sampling will contain the true parameter. As we can see in the Table 3 a 90% credible interval is given for each parameter of the model.
5.2. The estimated results figure
The sported plot of the estimated values are given as Figure 3;
77 79 SI S3 85 0 1000 2000 3000 4000
Figure 3: Displays the supported of the estimated parameters. The probability density functions of the marginalised posterior distributions, which are
calculated from the draws after warm-up in each chain, are shown on the left column side. The right column, which additionally displays the trace plots following warm-up. These plots show that all of the chains are oscillating around a set of values after the warm-up, which indicates that the MCMC has converged.
6. Model Checking
A common procedure in Bayesian statistics for validating your model is the posterior predective check [22]. We must first talk about the posterior predictive distribution in order to comprehend the posterior predictive check. The posterior predictive distribution in a mathematical notation.
v(y\y) = f y)v(6\y)d6 (21)
j 6
From the Equation 21, 6 is the parameter, in our posterior distribution which is discussed in the Section 4.4, 6 is the set of all parameters 6=(e, a, 61). The formula is merely provided for your convenience because we won't be using it to make posterior predictions. As an alternative, we will employ simulations to roughly approximate the posterior predictive distribution of y, which is really an extension of the posterior distribution of 6's approximation. The posterior predictive check simply makes a comparison between the observed data and the model's prediction we proceed to review the the graphical posterior predictive distribution of the total fertility rate.
Posterior Predict
1990 2000 2010 2020
time
Figure 4: Displays the posterior predicted.
From Figure 4, we can see that; the fitted model on the generated draws y — hat is very good fit. Therefore, we can results, the fitted model is a proper model for the given data-set.
6.1. Check Residuals
We may examine the residuals' autocorrelation using "checkresiduals()" in the forecast package will automatically run a test and display a few common diagnosis charts. Figure 5 is the common diagnostics charts of the total fertility rate.
Expected Values of the Posterior Predictive Errors
Figure 5: expected values of the posterior predictive errors.
The residuals series (upper part) shows that the given series is white noise. The dash blue lines are indicated the correlations are significantly not zero. In a white noise series, we expect each correlation to be zero. For a white noise series, we expect 95% of the spikes to lie within the where T is the length of the series. If more than 5% of the spikes are not in the blue bounds, then the series is not white noise probably. The histogram and quantile graph (middle part) show that the model is around normally distributed. Finally, the residual auto-correlation presents have showed that there is no correlation between the series segments. A summary result from the above interpretation is; we can step ahead and do our prediction based, as all the evidences present a good and fair series.
7. Prediction
Forecasting can be used in many situations: for making a schedule to the staff of a call center for some times in the future, we need to do a forecast as per call going to be received. For an inventory or stocking there is need for forecast of stock requirement. Forecast can be for few minutes beforehand till several years in advance. For making an efficient and effective plan, need to do forecasting. Forecasting can be done very easy and sometimes very tough. predictability for an event or quantity related to several factors;
• Understanding the factors to contribute on the event.
• How much data are accessible.
• Is there any effect of forecasting to the things are going to be forecasted [23].
Forecasting fertility; as the demographic characteristics (fertility, mortality, migration) of the human population is significant for socio-economics planning. There are mentioned lots of fertility forecasting methods, the methods discussed on the Section 1.1. Mostly, they have developed for the countries which are rich today and already passed their fertility transition, and at the current time they have low fertility patterns. Another point, is the developed methods are applied from the classical perspectives. A Bayesian projection model for the total fertility rate of country-specified projections of all countries have been done, the developed model consisted of three phases for the changes in total fertility rate, pre-transition (high fertility), the fertility
transition, and post-transition (low fertility). The model has been built on the United Nations Population Division's current [24].
In the same way, we aim to model the Total Fertility Rate of Afghanistan via the ARIMA model and achieve it is 16 years forecast, by harnessing the attributes of the Bayesian paradigm.
7.1. Forecasted Results
In the Table 4, which is numerical forecasted years, from 2022-2037 the 16 years of the total fertility rate of Afghanistan.
Table 4: Future forecasting
Point Forecast Lo.0.9 Hi.0.9 Lo.0.95 Hi.0.95
2023 3.853670 3.802970 3.903317 3.794032 3.915283
2024 3.718556 3.668241 3.768550 3.656223 3.776726
2025 3.586596 3.533241 3.639957 3.520394 3.650888
2026 3.455740 3.396801 3.510490 3.386231 3.520109
2027 3.332134 3.261010 3.388971 3.242005 3.400166
2028 3.207198 3.120882 3.278480 3.100771 3.290015
2029 3.085855 2.983132 3.162357 2.957278 3.173630
2030 2.967033 2.836914 3.055491 2.812259 3.068283
2031 2.851968 2.693900 2.959655 2.657657 2.972596
2032 2.739253 2.551682 2.859017 2.500849 2.878114
2033 2.626087 2.400823 2.762809 2.352485 2.776872
2034 2.519611 2.262385 2.677497 2.209586 2.693992
2035 2.411963 2.088992 2.588108 2.044437 2.604896
2036 2.311417 1.955690 2.510772 1.883435 2.533944
2037 2.210127 1.821591 2.432696 1.725132 2.452562
2038 2.114983 1.655969 2.371578 1.576980 2.388008
The variable's values we are are interested to forecast is unknown, therefore, it is a random variable. The random variable can take a range of random values. In Table 4, the first column is the years of forecast, the second column is the forecasted point or mean of the future values of the total fertility rate (which is the random variable). From columns 3-6 are the credible interval of 90% to 95%. For more details about the credible interval Section 5.1. The value for the total fertility rate of Afghanistan which is under consideration is given in Table 1, form 1982-2022. Afterwards, the value of the fertility of the country is random variable, because it is unknown and we are going to do forecast. In a Bayesian prediction we reach to the forecasting point which is the Table 4 is the numerical values of the forecasted amount.
7.2. Forecasted Results Figure
From Table 4, which is a 16 years Bayesian prediction of the total fertility rate of Afghanistan. Year of 2038 is showing a very good results. 2.1 is the Replacement fertility level: Total fertility of approximately 2.1 children per woman. This value shows the average number of children a woman would need to have in order to reproduce herself by having a daughter who reaches reproductive age. Now it is time to go to the graphical representation of the results. Figure 6, is going ot support the numerical values of Table 4, which is the forecasted values of the total fertility rate of Afghanistan from 2023-2038.
Figure 6: Future prediction
Figure 6 are consisted of two parts, from 1982-2022 the observed data, and 2023-2038 is the forecasted values of the total fertility rate of Afghanistan. The blue line in the Figure 6 shows the forecasted part of the value of the total fertility rate, which is the averages possible future value, and numerically was shown in the Table 4. In the Figure 6 it is obvious, there is two types of shadowed area, the 90% and the 95% credible intervals of the forecasted values. Forecast accuracy also can be analyzed form the shadowed areas, as much as, our forecasted time increase the credible interval or probability interval increase. And on the same way as much as years of forecast is near to the observed value the credible interval is more narrow and it shows the certainty or accuracy of the forecast.
8. Conclusion
In this article, we were started the work for doing Bayesian data analysis, using the total fertility rate of Afghanistan data-set. We discussed selecting a statistical model that is assumed to have generated the observed data, defining the prior and likelihood for the data, summarizing the posterior distribution, and using the posterior predictive check to assess the model's fit to the observed data. The model For the model ARIMA (1,2,1) has been chosen to the data-set. Model checking conducted, and 16 years forecasting has been done to the data. The total fertility rate of the country is projected to decline, and specifically it has been expected that approximately 2.1 (average number of children per woman) for the year 2038, this results is a very fundamental results that is the normal and standard measure of the fertility rate which is called replacement level of fertility. The importance of the demographic changes, obviously is a concern for the countries, for the effects on all affairs of human activity, economically, culturally, socially, and politically, therefore; the predicted value of the total fertility could be a good source for the government of the country for taken effective steps to control and decrease the fertility rate.
Need of the study: Afghanistan is a deve1oping country The country has been suffering for a long time, for a long-term war perspective, awareness of this type of issue could be significant for the government and policymakers, and non-government organizations to plan in order to monitor the fertility rate of the country. There was no such study, so there is a great need for the study, it
can be a good hint for those, who are going to work on the demographic area of Afghanistan.
One important tip is needed to mention is that we have done a prediction for the total fertility rate of the country for a specific period of future time, based on the series period of past time, no prediction can be a real future, but it is a motivation work for those, who are concern in such study.
References
Alkema, L., Raftery, A. E., Gerland, P., Clark, S. J., Pelletier, F., Buettner, T., and Heilig, G. K. (2011). Probabilistic projections of the total fertility rate for all countries. Demography, 48(3):815-839.
Bolanle, A. (2017). Forensic accounting as a tool for fighting corruption and terrorism: Case study of nigeria. Int. J. Soc. Sci. Econ. Res, 2(1):2210-2225.
Booth, H. (2006). Demographic forecasting: 1980 to 2005 in review. International journal of
forecasting, 22(3):547-581.
Box, G. E., Jenkins, G. M., Reinsel, G. C., and Ljung, G. M. (2015). Time series analysis: forecasting and control. John Wiley, Sons.
Carmichael, G. A. et al. (2016). Fundamentals of demographic analysis: Concepts, measures, and methods. Springer.
Gamerman, D. and Lopes, H. F. (2006). Markov chain Monte Carlo: stochastic simulation for Bayesian inference. CRC press.
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2013). Bayesian data analysis. CRC press.
Gelman, A., Hwang, J., and Vehtari, A. (2014). Understanding predictive information criteria for bayesian models. Statistics and computing, 24(6):997-1016.
Gelman, A., Meng, X.-L., and Stern, H. (1996). Posterior predictive assessment of model fitness via realized discrepancies. Statistica sinica, pages 733-760.
Hoffman, M. D., Gelman, A., et al. (2014). The no-u-turn sampler: adaptively setting path lengths in hamiltonian monte carlo. J. Mach. Learn. Res., 15(1):1593-1623. Holmes, E. E., Scheuerell, M. D., and Ward, E. (2020). Applied time series analysis for fisheries and environmental data. NOAA Fisheries, Northwest Fisheries Science Center, Seattle, WA.
Hyndman, R. J. and Athanasopoulos, G. (2018). Forecasting: principles and practice. OTexts. Knudsen, C., McNown, R., and Rogers, A. (1993). Forecasting fertility: An application of time series methods to parameterized model schedules. Social Science Research, 22(1):1-23. Korner-Nievergelt, F., Roth, T., Von Felten, S., Gu?elat, J., Almasi, B., and Korner-Nievergelt, P. (2015). Bayesian data analysis in ecology using linear models with R, BUGS, and Stan. Academic Press.
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. Lee, R. D. and Carter, L. R. (1992). Modeling and forecasting us mortality. Journal of the American statistical association, 87(419):659-671.
Lee, R. D. and Tuljapurkar, S. (1994). Stochastic population forecasts for the united states: Beyond high, medium, and low. Journal of the American Statistical Association, 89(428):1175-1189.
Marriott, J., Ravishanker, N., and Gelfand, A. E. (1993). Bayesian analysis of arma processes: Complete sampling based inference under full likelihoods. Technical report, STANFORD UNIVCA DEPT OF STATISTICS.
McDonald, J. (1979). A time series approach to forecasting australian total live births. Demography, 16(4):575-601.
McNown, R. and Rogers, A. (1989). Forecasting mortality: A parameterized time series approach.Demography, 26(4):645-660.
[21] Ortega, J. A. and Poncela, P. (2005). Joint forecasts of southern european fertility rates with non-stationary dynamic factor models. International Journal of Forecasting, 21(3):539-550.
[22] Raziq, A., Sheikh, N., Rafique, M., Murtaza, G., and Muhammad, A. (2021). Forecasting fertility and child mortality rates of pakistan using time series model. Linguistica Antverpiensia, 2021(2):1194-1213.
[23] Tripathi, P. K., Mishra, R. K., and Upadhyay, S. (2018). Bayes and classical prediction of total fertility rate of india using autoregressive integrated moving average model. Journal of Statistics Applications , Probability, 7(2):233-244.
[24] Watanabe, S. and Opper, M. (2010). Asymptotic equivalence of bayes cross validation and widely applicable information criterion in singular learning theory. Journal of machine learning research, 11(12).
