Next Article in Journal
Dynamics and the Cohomology of Measured Laminations
Previous Article in Journal
Cost Effectiveness Analysis of Optimal Malaria Control Strategies in Kenya
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

New Method of Randomized Forecasting Using Entropy-Robust Estimation: Application to the World Population Prediction

by
Yuri S. Popkov
1,2,3,*,
Yuri A. Dubnov
1,2 and
Alexey Yu. Popkov
1,2
1
Institute for Systems Analysis of Russian Academy of Sciences, Moscow 117312, Russia
2
Moscow Institute of Physics and Technology (State University), Moscow 141700, Russia
3
National Research University Higher School of Economics, Moscow 101000, Russia
*
Author to whom correspondence should be addressed.
Mathematics 2016, 4(1), 16; https://doi.org/10.3390/math4010016
Submission received: 13 October 2015 / Revised: 25 January 2016 / Accepted: 4 March 2016 / Published: 11 March 2016

Abstract

:
We propose a new method of randomized forecasting (RF-method), which operates with models described by systems of linear ordinary differential equations with random parameters. The RF-method is based on entropy-robust estimation of the probability density functions (PDFs) of model parameters and measurement noises. The entropy-optimal estimator uses a limited amount of data. The method of randomized forecasting is applied to World population prediction. Ensembles of entropy-optimal prognostic trajectories of World population and their probability characteristics are generated. We show potential preferences of the proposed method in comparison with existing methods.

1. Introduction

For a studied process, forecasting as a procedure consists of four consecutive stages: modeling (model design), learning (estimation of model characteristics), testing (of the “learned” model) and prediction of future development.
Forecasting is based on retrospective data analysis with its subsequent extrapolation to future periods. Consider the state of a studied process at moment t 0 , and suppose that the problem is to forecast further evolution of the process on a time interval T f r c = [ t 0 , T ] . Then, it is necessary to operate existing data on its past dynamics on a time interval T r t s = [ T , t 0 ] , where T < t 0 (the so-called retrospective data). Generally, retrospective data and the time interval T r t s are divided into two groups, namely data serving for the estimation of the model’s characteristics (on a time interval T e s t = [ T , t e ] ) and model testing (on a time interval T t s t = [ t e , t 0 ] ).
There exist at least three forecasting techniques differing in the objectification degree of constructed forecasts. The first technique, referred to as scenario forecasting [1], proceeds from the scenario approach whose objectification is replaced by the opinion of an expert group. Actually, it implements only the stages of modeling and prediction: learning and testing are eliminated owing to the opinion of experts, who choose an appropriate mathematical model of a studied process and form value sets (scenarios) of the model parameters. Then, the model with the scenario parameter values generates forecasting trajectories. As a matter of fact, such a forecasting technique is most widespread in demographic prediction [2,3,4]. Note that real retrospective data about a studied process are not utilized. They are indirectly reflected by the knowledge and experience of invited experts.
The second technique of forecasting explicitly involves real data for model learning and testing. The framework of mathematical statistics provides numerous estimation methods for model parameters; see [5,6,7,8]. Here, a major assumption is that the model possesses deterministic parameters, the values of which are defined using sets of real retrospective data. The latter are treated as a stochastic object with certain properties (a sample from a universe, normal distribution, etc.). In this case, one may assign different probabilistic characteristics (variances, confidence intervals, and so on) to the derived estimates of model parameters.
Sometimes, these characteristics assist in constructing the probabilistic characteristics of forecasting trajectories. The described technique will be termed probabilistic forecasting (PF) [9,10,11,12,13,14]. We emphasize that the above hypotheses regarding the stochastic properties of real datasets are almost impossible to verify, especially under small data arrays. It follows for the low efficiency of forecasts [15,16].
Finally, the third technique of forecasting proposed in this paper stems from the randomized model (RM) of a studied process, where model parameters are supposed random. Hence, we characterize RMs using the probability density functions (PDFs) of their parameters. At the learning stage, the PDFs of the model parameters are estimated on the basis of real retrospective data.
A randomized model generates an ensemble of forecasting trajectories, where each trajectory corresponds to a set of random realizations of parameters with the derived estimates of the PDFs. Computer simulation of such models employs the Monte Carlo method. Below, this technique will be called randomized forecasting (RF).
As opposed to existing methods [17,18,19], the proposed method of randomized forecasting is based on entropy-optimal estimations of PDFs for real datasets. The structure of the randomized dynamic model, used in this method, is based on ordinary differential equations.
The developed method serves for obtaining randomized predictions of the World population dynamics. Modeling the World population variations in time and space forms a major problem of demographic analysis [20,21,22].
Throughout the paper, we describe the above mentioned dynamics by the exponential model incorporating several parameters associated with fertility and mortality rates, as well as its change in time. In the randomized setting, they are assumed random, whereas World population is measured with random errors. To find the corresponding PDFs, the method involves the retrospective population data provided by the UN (see UNdata service at https://data.un.org/). In addition, we perform the comparative analysis of the PF- and RF-based approaches.

2. Randomized Model: Linear Differential Form

Consider a dynamic object having an input f ( t ) = { f 1 ( t ) , , f m ( t ) } and an output x ( t ) = { x 1 ( t ) , , x n ( t ) } . The components of the input and output vectors can be observed (measured) on a time interval T r t s = [ T , t 0 ] .
The relationship between the input and output of the object is described by the linear nonautonomous system of ordinary differential equations:
d x ( t ) d t = A x ( t ) + B f ( t ) , x ( T ) = x 0
where x R n ; f R m , m n ; A = [ a i j | ( i , j ) = 1 , n ¯ ] and B = [ b i k | i = 1 , n ¯ , k = 1 , m ¯ ] denote matrices of appropriate dimensions.
The object’s output is observed with inevitable disturbances modeled by a vector noise ξ ¯ ( t ) = { ξ 1 ( t ) , , ξ n ( t ) } . Therefore, the observed output of the model acquires the form:
v ( t ) = x ( t ) + ξ ¯ ( t )
where v ( t ) R n .
Equations (1) and (2) define a linear dynamical randomized model (LDRM) if:
  • the matrix A is a random matrix (with independent random elements or elements with independent random components) of the interval type:
    A = [ A : A A A + ]
    where A , A + mean given matrices;
  • there exists a probability density function (PDF) P ( A ) , A A ;
  • the vector ξ ¯ is a random vector (i.e., contains independent random components) of the interval type:
    Ξ = [ ξ ¯ ξ ¯ ξ ¯ + ]
  • there exists a probability density function (PDF) Q ( ξ ¯ ) , ξ ¯ Ξ ;
  • the matrix B possesses fixed known elements.
Under the stated conditions, the LDRM generates an ensemble of random trajectories on the time interval T = [ T , T ] .
Let us rewrite the LDRM (2) in the input-output representation using the notion of matrix exponent [23]:
W ( A | t τ ) = exp [ A ( t τ ) ]
The input and output are measured at discrete moments with step h. Hence, on the time interval T e s t , we have:
x [ T + i h ] = W ( A | i h ) x 0 + + T T + i h W ( A | T + i h , τ ) B f ( τ ) d τ , i 0 , N e s t ¯
here N e s t = [ ( t e T ) / h ] and [ ] indicates the integer part of •.
The LDRM output Equation (3) observed at discrete moments has the following form:
v [ T + i h ] = x [ T + i h ] + ξ ¯ [ T + i h ] , i = 0 , N e s t ¯
Let us denote ξ ¯ ( i ) = ξ ¯ [ T + i h ] , i = 0 , N e s t ¯ . They are random vectors with independent and interval components. Now, we introduce the block-vector ξ ^ = { ξ ¯ ( 0 ) , , ξ ¯ ( N e s t ) } . As we assume that these vectors and their components are independent, then the joint PDF is:
Q ( ξ ^ ) = i = 0 N e s t Q i ( ξ ¯ ( i ) ) = i = 0 N e s t j = 1 n q i j ( ξ j ( i ) )
The domain of this function is:
Ξ ^ = Ξ × Ξ × Ξ ( N e s t + 1 ) multipliers
Here, ξ ¯ ( T ) , ξ ¯ ( T + h ) , , ξ ¯ ( T + N e s t h ) gives a sequence of n-dimensional independent random vectors of the interval type, associated with corresponding PDFs.
As soon as the matrix A and the noise vector ξ ^ are random and characterized by the PDFs P ( A ) and Q ( ξ ¯ ) , respectively, so the observed output of the LDRM represents an ensemble V of random trajectories v [ T + i h ] , i = 0 , N e s t ¯ .

3. S PQ 1 Entropy-Robust Estimation

The first stage of the randomized forecasting (RF) is an estimation of the PDFs of the RMs’ parameters and of the measurement noises. It is a classical problem of the Bayesian approach, and there exists the classical maximum likelihood method (or maximum relation of likelihood; see [6]) for its solving.
Let us recall the definition of the function of the relation of likelihood (FRL) in the terms of Section 2. The a priori PDFs P 0 ( A ) , Q 0 ( ξ ^ ) and the a posteriori PDFs P ( A ) , Q ( ξ ^ ) are the basic notations of the FRL. The FRL takes the form:
L ( A ) = ln P ( A ) P 0 ( A ) , L ( ξ ^ ) = ln Q ( ξ ^ ) Q 0 ( ξ ^ )
If the PDFs in these expressions are restored as functions of the parameters, then maximization of these functions gives “optimal” estimations. The principle of maximization of the FRL takes the form:
A ^ = arg max A L ( A )
As a declaration, this principle is fine. However, how is it possible to restore the PDFs P ( A ) and P 0 ( A ) (also Q ( ξ ^ ) and Q 0 ( ξ ^ ) )? The problem of the restoration of the PDFs remains outside of the FRL.
Let us consider the functional of the likelihood relation (FuRL) in the following form:
L [ P ( A ) ] = A P ( A ) ln P ( A ) P 0 ( A ) d A , L [ Q ( ξ ^ ) ] = Ξ ^ Q ( ξ ^ ) ln Q ( ξ ^ ) Q 0 ( ξ ^ ) d ξ ^
From Equation (12), we can see that the FuRL is the mathematical expectation of the FRL. On the other side, the FuRL is the opposite generalized information Boltzmann entropy (Kullback–Leibler distance) [24,25], that is:
H A [ P ( A ) ] = L [ P ( A ) ] , H [ Q ( ξ ^ ) ] = L [ Q ( ξ ^ ) ]
According to [26], maximization of entropy functions gives the best robust solution under high uncertainty. This idea with the addition of real data balance conditions forms the basis of the S P Q 1 entropy-robust estimation method [27].
The S P Q 1 entropy-robust estimation can be reformulated as a problem of functional nonlinear programming [28,29]:
H [ P ( A ) , Q ( ξ ¯ ) ] = A P ( A ) ln P ( A ) d A + + Ξ ¯ Q ( ξ ^ ) ln Q ( ξ ^ ) d ξ ^ min
subject to the constraints imposed on:
-
the class of (normalized) PDFs:
A P ( A ) d A = 1 Ξ ^ Q ( ξ ^ ) d ξ ^ = 1
-
the balance between the first moment vector of the observed output v [ T + i h ] = v ( i ) in the LDRM Equation (7) and the real data vector y [ T + i h ] = y ( i ) :
M { v ( i ) } = v ¯ ( i ) [ P ( A ) , Q ( ξ ^ ) ] = A u ( i ) ( A ) P ( A ) d A + + Ξ ^ ξ ¯ ( i ) Q ( ξ ^ ) d ξ ^ = y ( i ) , i 0 , N e s t ¯
where:
u ( i ) ( A ) = W ( A | i h ) x 0 + T T + i h W ( A | T + i h τ ) B f ( τ ) d τ
S P Q 1 entropy-robust estimation uses the first moment vector of the output LDRM. It is possible to use moments of higher order. This depends on measurable real data. If they represent the k-moment, then the balance condition can be formulated in the following form:
M { v k } 1 / k = y
where v ( k ) is a vector of k-moment components. In this case, we will have S P Q k entropy-robust estimation.
The problem Equations (14)–(17) are related to the Lyapunov problem [28,29] where the goal functional and constraints are of an integral type. Here, we will use the necessary condition of equality to zero of the integral equation:
X h ( x ) g ( x ) d x = 0
where the function h ( x ) is continuous and is equal to zero on the boundary of the set X (the class C ˜ ); the function g ( x ) is differentiable (the class D ). Then, this equality will be valid for any function h ( x ) with the mentioned properties if:
g ( x ) = 0 , x X
Now, we return to the problem Equations (14)–(17) and introduce the Lagrange functional:
L [ P ( A ) , Q ( ξ ^ ) ] = H [ P ( A ) , Q ( ξ ^ ) ] + λ A P ( A ) d A 1 + μ Ξ ^ Q ( ξ ^ ) d ξ ^ 1 + + i = 0 N e s t θ ¯ ( i ) , m ( i ) [ P ( A ) , Q ( ξ ^ ) ]
where:
m ( i ) [ P ( A ) , Q ( ξ ^ ) ] = v ¯ ( i ) [ P ( A ) , Q ( ξ ^ ) ] y ( i )
Sign , denotes a scalar product.
As the solution of the problem Equations (14)–(17) is searched in the class of differentiable functions, then the Gato derivation can be used for determination of the variation of Lagrange functional Equation (21).
Let us denote the solution of the problem Equations (14)–(17) as P * ( A ) and Q * ( ξ ^ ) . Furthermore, introduce the functions ϕ ( A ) C ˜ , ψ ( ξ ^ ) D and two scalar variables α , β to present the functions P ( A ) , Q ( ξ ^ ) in the following form:
P ( A ) = P * ( A ) + α ϕ ( A ) , Q ( ξ ^ ) = Q * ( ξ ^ ) + β ψ ( ξ ^ )
Functions P * ( A ) , Q * ( ξ ^ ) , as a solutions of the problem Equations (14)–(17), are fixed. The optimality conditions for the problem Equations (14)–(17) take the form:
d L d α α = β = 0 = 0 , d L d β α = β = 0 = 0
The application of these conditions leads to the following systems of integral equations:
A ϕ ( A ) ln P ( A ) + 1 + λ + i = 0 N e s t θ ^ ( i ) , u ( i ) ( A ) d A = 0 Ξ ^ ψ ( ξ ^ ) ln Q ( ξ ^ ) + 1 + μ + i = 0 N e s t θ ¯ ( i ) , ξ ¯ ( i ) d ξ ^ = 0
According to Equations (19) and (20), we obtain the following equations, which are necessary optimality conditions (necessary conditions of Lagrangian-stationarity) for the problem Equations (14)–(17):
ln P ( A ) + 1 + λ + i = 0 N e s t θ ^ ( i ) , u ( i ) ( A ) = 0 ln Q ( ξ ^ ) + 1 + μ + i = 0 N e s t θ ¯ ( i ) , ξ ¯ ( i ) = 0
The solution of the problem Equations (14)–(17) takes the form:
P * ( A ) = exp i = 0 N e s t θ ¯ ( i ) , u ( i ) ( A ) R θ ¯ ( 0 ) , , θ ¯ ( N e s t ) Q * ( ξ ^ ) = exp i = 0 N e s t θ ¯ ( i ) , ξ ¯ ( i ) Q θ ¯ ( 0 ) , , θ ¯ ( N e s t )
where:
R θ ¯ ( 0 ) , , θ ¯ ( N e s t ) = A exp i = 0 N e s t θ ¯ ( i ) , u ( i ) ( A ) d A Q θ ¯ ( 0 ) , , θ ¯ ( N e s t ) = Ξ ^ exp i = 0 N e s t θ ¯ ( i ) , ξ ¯ ( i ) d ξ ¯ ( 0 ) d ξ ¯ ( N e s t )
The vectors of Lagrange multipliers are determined from the following equations:
U ( θ ¯ ( 0 ) , , θ ¯ ( N e s t ) ) = A u ( i ) ( A ) exp i = 0 N e s t θ ¯ ( i ) , u ( i ) ( A ) R θ ¯ ( 0 ) , , θ ¯ ( N e s t ) d A + + Ξ ^ ξ ¯ ( i ) exp i = 0 N e s t θ ¯ ( i ) , ξ ¯ ( i ) Q θ ¯ ( 0 ) , , θ ¯ ( N e s t ) d ξ ¯ ( 0 ) d ξ ¯ ( N e s t ) y ( i ) = 0 , i [ 0 , N e s t ]
Calculation of the vectors θ ^ * = { θ ¯ * ( 0 ) , , θ ¯ * ( N e s t ) } is turned into a search for the global minimum of the residual function:
J ( θ ^ ) = U ( θ ^ ) L 2
A global optimization algorithm is based on the simple Monte Carlo trials proposed in [30]. However, as soon as the L 2 metric is a convex function, one of the traditional gradient-based local optimization methods can be used for its solving.

4. Randomized Forecast Implementation

We comprehend a randomized forecast as an ensemble of trajectories on a forecasting interval T f r c = [ t 0 , T ] , which has to be generated using the model Equations (6) and (7) with the random matrix A and noise ξ ¯ described by the PDFs P * ( A ) and Q * ( ξ ¯ ) , respectively; see formulas Equations (27) and (28). The matrices A and the vector noise ξ ¯ belong to the parallelepipeds A from Equation (3) and Ξ from Equation (4), respectively.
Let us study the generation problem of random matrices with the PDF Equation (27). First, we transform a matrix into a vector through concatenation of its rows. This procedure yields a vector a of length m = n 2 . Additionally, the domain of random matrices becomes an m-dimensional parallelepiped:
A = [ a a a + ]
where the vectors a and a + result from the row concatenation of the matrices A and A + , respectively.
Lets consider a transformation of a vector a into a vector q belonging to the m-dimensional unit nonnegative cube Q :
a = ( a + a ) q + a , Q = { q : 0 q 1 }
Therefore, the entropy-optimal PDF undergoes the following chain of transformations:
P * ( A ) P ( a ) P ( q )
Therefore, it is necessary to generate random vectors q Q with PDFs P ( q ) . The generation was implemented by the acceptance-rejection algorithm [24].

5. Application of the RF Method for World Population Prediction

5.1. The World Population Prediction Problem

The state of an isolated population is characterized by its size E ( t ) on a calendar time interval T = [ T , t 0 ] . Population size varies under the impact of fertility and mortality processes, since World population is an isolated system. Within the framework of the linear population dynamics model, fertility and mortality are described by corresponding rates, whereas the flows of newborns and decedents appear proportional to population size, while fertility (b) and mortality (m) rates are considered as linear time-dependent parameters [22].

5.1.1. Randomized population model

World population evolves according to the following differential equation that has an analytical solution:
d E ( t ) d t = ( b ( t ) m ( t ) ) E ( t ) , b ( t ) = b 0 + u b t , m ( t ) = m 0 + u m t E ( t ) = E 0 exp ( b 0 m 0 ) + 1 2 ( u b u m ) t t , where E 0 = E ( T ) , b 0 = b ( T ) , m 0 = m ( T )
Real measurements of the population size dynamics modeled by Equation (34) take place at discrete moments. Hence, the population size at discrete moments i h (where h specifies a given increment) is defined by the expression:
Φ i ( r , u r | E 0 ) = E 0 exp ( r + u r i ) i h , i I r = b 0 m 0 , u r = 1 2 u b u m h
with parameters r and u r , which describe the result of the difference between fertility and mortality flows.
World population is measured in billions of people. Fertility and mortality processes aggregate many factors, including measurement errors, whose quantitative analysis is impossible or complicated. On the other hand, the mass nature of fertility and mortality processes admits their modeling based on the probabilistic approach.
Thus, the resultant flow rate and its changing in time are supposed to be random variables with a joint probability density function P ( r , u r ) defined on the rectangle:
J = I r I u r , I r = [ r , r + ] , I u r = [ u r , u r + ]
Generally, measurement errors are modeled by an additive noise ξ [ i h ] of the interval type:
ξ [ i h ] Ξ i = [ ξ i , ξ i + ] , i I
where I indicates a time interval of such measurements. By assumption, the PDFs q i ( ξ [ i h ] ) , i I are specified on the intervals Ξ i from Equation (37). Due to the independence of the set of random variables ξ [ 0 ] , , ξ [ i h ] , their joint PDF has the form:
Q ( ξ ¯ ) = i I q i ( ξ [ i h ] )
Therefore, the randomized model of World population dynamics can be described by:
v [ i h ] = Φ i ( r , u r | E 0 ) + ξ [ i h ] , i I
where the function Φ i ( r , u r | E 0 ) meets Equality (35).

5.1.2. Real and forecasting data

For PDF estimation, address the World population measurements for the period from 1960–1995 (http://data.un.org/; see Table 1).
The entropy-optimal RM is tested via the measurements of the World population dynamics during the period from 1995–2015 and the values E 1985 p r n for this period (according to the UN forecast announced in 1985) (http://www.irbis.vegu.ru/repos/1002/Html/27.htm; see Table 2).
The World population prediction till 2050 using the UN data is illustrated by Table 3 (http://data.un.org/; [31]). United Nations’ projections both for the testing interval and the forecasting interval were made in accordance with the commonly-used cohort-component method, based on age-specific estimates of the components of population change (fertility, mortality and international migration) [32,33]. We will compare this prediction with its randomized analog.
World population is measured in billions of people. In the sequel, the subscript r e a l indicates the measured data of the population.
To summarize, the RM has the following forms on corresponding time intervals:
  • on the estimation interval T e s t (see Table 1):
    v [ i h ] = Φ i ( r , u r | E r e a l e s t [ 0 ] ) + ξ [ i h ] , i [ 0 , 7 ]
  • on the testing interval T t s t (see Table 2):
    v [ i h ] = Φ i ( r , u r | E r e a l t s t [ 0 ] ) + ξ [ i h ] , i [ 0 , 4 ]
  • on the forecasting interval T f r c (see Table 3):
    v [ i h ] = Φ i ( r , u r | E r e a l f r c [ 0 ] ) , i [ 0 , 5 ]
where r and u r are random parameters with the entropy-optimal PDFs P * ( r , u r ) and ξ ¯ is a vector of random noise with entropy-optimal PDF Q * ( ξ ¯ ) .

5.1.3. The entropy-optimal PDFs of the parameters and noise

According to the general entropy-robust estimation procedure of PDFs (see Section 2), we have:
  • the PDF of the parameters r and u r in the form:
    P * ( r , u r ) = 1 R ( θ ¯ | E r e a l e s t [ 0 ] ) j = 0 7 p j * ( r , u r | θ j ) , p j * ( r , u r | θ j ) = exp θ j Φ j ( r , u r | E r e a l e s t [ 0 ] )
  • the PDF of the noise in the form:
    Q * ( ξ ¯ ) = 1 Q ( θ ¯ ) j = 0 7 q j * ( ξ [ j h ] | θ j ) , q j * ( ξ [ j h ] | θ j ) = exp θ j ξ [ j h ]
    where:
    R ( θ ¯ | E r e a l e s t [ 0 ] ) = J j = 0 7 exp θ j Φ j ( r , u r | E r e a l e s t [ 0 ] ) d r d u r
    Q ( θ ¯ ) = j = 0 6 ξ j ξ j + exp ( θ j ξ [ j h ] ) d ξ [ j h ] = = j = 0 6 1 θ j exp ( θ j ξ j ) exp ( θ j ξ j + )
To calculate the Lagrange multipliers, we will solve the system of balance equations (see Equations (16) and (17):
1 R ( θ ¯ | E r e a l e s t [ 0 ] ) J Φ i ( r , u r | E r e a l e s t [ 0 ] ) j = 0 7 exp θ j Φ j ( r , u r | E r e a l e s t [ 0 ] ) d r d u r + + 1 Q ( θ ¯ ) Ξ ξ [ i h ] j = 0 7 exp θ j ξ [ j h ] d ξ [ i h ] E r e a l e s t [ i h ] = 0 , i [ 0 , 7 ]
We will denote:
N i ( θ ¯ | E r e a l e s t [ 0 ] ) = J Φ i ( r , u r | E r e a l e s t [ 0 ] ) j = 0 7 exp θ j Φ j ( r , u r | E r e a l e s t [ 0 ] ) d r d u r
Then, Equation (47) can be rewritten as:
G i ( θ ¯ | E r e a l e s t [ 0 ] ) = N i ( θ ¯ | E r e a l e s t [ 0 ] ) R ( θ ¯ | E r e a l e s t [ 0 ] ) + L i ( θ i ) y i = 0 y i = E r e a l e s t [ i h ] , i = [ 0 , 7 ]
where:
L i ( θ i ) = exp ( θ i ξ i ) ( ξ i + 1 θ i ) exp ( θ i ξ i + ) ( ξ i + + 1 θ i ) exp ( θ i ξ i ) exp ( θ i ξ i + )
We endeavor to solve these equations through minimizing the residual function:
J ( θ ¯ ) = G ( θ ¯ ) min

5.2. The Results of Computer Experiment

5.2.1. General conditions

We have performed calculations for the stages of estimation, testing and prediction with the following ranges for the model parameters, which include the possibility for both positive and negative trends of World population growth:
I r = [ 0 . 025 ; 0 . 075 ] , I u r = [ 0 . 002 ; 0 . 001 ]
and identically the range for the measurement noises:
Ξ j = [ 0 . 5 ; 0 . 5 ] , j [ 0 , 7 ]
Under the above ranges of the model parameters and measurement noises, the computer experiment has constructed the entropy-optimal PDFs and generated the ensembles of corresponding random trajectories for RM testing and randomized forecasting.

5.2.2. Estimation

On the estimation interval, we employ the data from Table 1. The residual function J ( θ ¯ ) Equation (51) is a function of eight variables, and it contains integral component Equations (45)–(49) to be estimated only numerically. For this, we have selected the so-called tiled method of two-dimensional integral estimation, which represents a combination of several quadrature formulas [34].
The idea of the tiled method consists of: (1) partitioning the whole domain of integration into a set of smaller-area subdomains having the rectangular or trapezoidal shape; and (2) applying appropriate quadrature formulas on each subdomain. The described method is implemented in MATLAB by the function quad2d.
Minimization of the residual function J ( θ ¯ ) (51) runs by the nonlinear trusted region method implemented by the function lsqnonlin from the package Optimization Toolbox. The function lsqnonlin has been optimized [35] for nonlinear least-squares problems.
The function lsqnonlin have several user-defined options for stopping criteria, such as function evaluation tolerance ( 10 6 ), step size tolerance ( 10 6 ) and maximum number of iterations (500). Table 4 presents the calculated Lagrange multipliers for the above ranges.
Figure 1 and Figure 2 demonstrate the entropy-optimal PDFs for the parameters of Model (35) and noise components. The functions P * ( r , u r , θ ¯ ) and q * ( ξ [ i h ] ) , i [ 0 , 7 ] represent an entropy-optimal distribution of random variables at the corresponding intervals and will be used for making randomized predictions.

5.2.3. Testing

Testing of the model has been made using the data from Table 2. The population size has been evaluated by the Formula (41), where r , u r are the random variables with the PDF P * ( r , u r ) and ξ [ i h ] are the random noises with the PDFs q i ( ξ [ i h ] ) , i [ 0 , 4 ] (see Figure 1 and Figure 2).
To generate an ensemble of random variables, the two-dimensional modification of the Ulam–von Neumann method has been used [24]. The size of the generated sample is k = 100,000.
Each pair of the random values r and u r defines a separate exponential growth curve; moreover, for each point i h , a random value of the noise ξ [ i h ] is added according to its PDF. As a result, the constructed trajectory of World population dynamics is not an exponential function.
The test procedure yields the probabilistic characteristics of the parameters r, u r and ξ. It depends on the initial population size for the testing interval. Figure 3 shows the ensemble of the model-based trajectories for the parameter range from Equation (52) and the noise range from Equation (53), with a data-based selection of the initial point:
E t s t [ 0 ] = E r e a l t s t [ 0 ]
Figure 3 has the following notation: 1, the ensemble-average trajectories of population dynamics; 2, real population dynamics on the testing interval; 3, population dynamics on the testing interval according to the UN forecast of the year 1998; 4, the boundaries of the first and the third quartiles. The same model-generated ensemble for testing interval can be presented as a box plot with median mark and interquartile ranges; see Figure 4.
Testing quality is assessed by the root-mean-square error of the average trajectory with respect to its real counterpart:
δ = E t s t E r e a l t s t = i = 0 4 ( E r e a l t s t [ i h ] E t s t [ i h ] ) 2
as well as by the relative error:
ε = δ E r e a l t s t + E t s t
For instance, the UN forecast for the testing interval (Table 2) has the following errors:
δ 1985 = 0 . 228 , ε 1985 = 0 . 008
In our case, the deviation of the model-average trajectory from the real one (Figure 3) appears to be appreciably smaller and demonstrates the following errors:
δ R M = 0 . 079 , ε R M = 0 . 003

5.2.4. Prediction

This simple RM has been applied to predict World population dynamics for the period from 2015–2050. The trajectory ensemble corresponding to the UN predictions (Table 3) is illustrated by Figure 5: 1, the ensemble-average trajectory; 2, UN projection; 4, the boundaries of the interquartile range (IQR zone). The corresponding box-plot is presented in Figure 6.
The presented results testify that the randomized forecasting as opposed to existent methods gives a set of probability characteristics of the World population prediction, which is calculated by using the ensemble of prognostic trajectories. The latter is generated by the randomized dynamic model with entropy-optimal PDFs of parameters. The randomized projection algorithm shows significantly closer to real data numerical results for testing interval, as well as stable projection that is higher, but close to the modern UN forecast for the future. According to our randomized model, 2026 will be the first eight billion year.

6. Conclusions

In this paper, we have suggested a randomized forecasting method that operates dynamic models described by linear ordinary differential equations with random parameters. Entropy-robust estimation has been developed for the probability density functions (PDFs) of model parameters and noisy measurements based on entropy maximization. It has been shown that the above PDFs belong to the exponential class. The randomized forecasting technique has been applied to randomized prediction of the World population dynamics. It has been demonstrated that randomized forecasting gives a set of probability characteristics of the World population dynamics.

Acknowledgments

This work was supported by the Russian Foundation for Basic Research (Project No. 16-07-00743). We also thank the comments from the three anonymous reviewers, which improved the quality of the paper.

Author Contributions

All authors have contribued significantly to the work reported by this article. Y.S. developed and declared the method; Y.S. and Y.A. conceived and designed the experiments; Y.A. performed the experiments; Y.A. and A.Y. analyzed the data; A.Y. contributed reagents/materials/analysis tools; Y.S., Y.A and A.Y wrote the paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. De Beer, J.; van Wissen, L. (Eds.) Europe: One Continent, Different Worlds. Population Scenarios for the 21th Century; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1999.
  2. Hilderink, H. World Population in Transition. An Integrated Regional Modelling Framework; Groningen University Press: Groningen, The Netherlands, 2000. [Google Scholar]
  3. Feichtinger, G. (Ed.) Vienna Yearbook of Population Reseach; Vienna Institute of Demography, Austrian Academy of Sciences: Vienna, Austria, 2003.
  4. Wissen, L.J.G. PROFILE: Population Research in the Netherlands. Public Serv. Rev. Eur. Union 2013, 25, 33. [Google Scholar]
  5. Kendall, M.G.; Stuart, A. The Advanced Theory of Statistics; Vol. II, Inference and Relationship, 3rd ed.; Griffin: London, UK, 1973. [Google Scholar]
  6. Cramer, H. Mathematical Methods of Statistics; Princeton Univ. Press: Princeton, NJ, USA, 1999. [Google Scholar]
  7. Aivazyan, S.A.; Mkhitryan, V.S. Prikladnaya statistika i osnovy ekonometriki. In Applied Statistics and Foundations of Econometrics; Interreklama: Moscow, Russia, 2003. [Google Scholar]
  8. Kobzar, A.I. Prikladnaya matematicheskaya statistika. In Applied Mathematical Statistics; Fizmatlit: Moscow, Russia, 2006. [Google Scholar]
  9. Cho, S. A Linear Bayesian Stochastic Approximation to Update Project Duration Estimation. Eur. J. Oper. Res. 2009, 196, 585–593. [Google Scholar] [CrossRef]
  10. Zellner, A. Bayesian Shrinkage Estimation and Forecasts of Individual and Total or Aggregate Outcomes. Econ. Model. 2010, 27, 1392–1397. [Google Scholar] [CrossRef]
  11. Horvath, R. Research & Development and Growth: A Bayesian Model Averaging Analysis. Econ. Model. 2011, 28, 2669–2673. [Google Scholar]
  12. Kim, M.J.; Jiang, R.; Makis, V.; Lee, C.-G. Optimal-Bayesian Fault Prediction Scheme for Partially Observable System Subject to Random Failure. Eur. J. Oper. Res. 2011, 214, 331–339. [Google Scholar] [CrossRef]
  13. Musal, R.M.; Soyer, R.; McCabe, C.; Kharroubi, S.A. Estimating the Population Utility Function: A Parametric Bayesian Approach. Eur. J. Oper. Res. 2012, 218, 538–547. [Google Scholar] [CrossRef]
  14. Borisov, A.V. A posteriori Minimax Estimation with Limit Likelihood. Autom. Remote Control 2012, 9, 49–56. [Google Scholar]
  15. Lawrence, M.; Goodwin, P.; O’Connor, M.; Öncal, D. Judgemental Forecasting: A Review of Progress over the Last 25 years. Int. J. Forecast. 2006, 22, 493–518. [Google Scholar] [CrossRef]
  16. Kociecki, A.; Kolasa, M.; Rubaszek, M. A Bayesian Method of Combining Judgmental and Model-Based Density Forecasts. Econ. Model. 2012, 29, 1349–1355. [Google Scholar] [CrossRef]
  17. Lahiri, K.; Peng, H.; Zhao, Y. Testing the Value of Probability Forecasts for Calibrated Combining. Int. J. Forecast. 2015, 31, 113–129. [Google Scholar] [CrossRef] [PubMed]
  18. Jacobson, M.Z. Fundamentals of Atmospheric Modeling; Cambridge University Press: New York, NY, USA, 2005; p. 828. [Google Scholar]
  19. Allen, M.R.; Stainforth, D.A. Towards Objective Probabilistic Climate Forecasting. Nature 2002, 419. [Google Scholar] [CrossRef] [PubMed]
  20. Popkov, A.Y.; Popkov, Y.S.; van Wissen, L. Positive Dynamic Systems with Entropy Operator: Application to Labour Market Modeling. Eur. J. Oper.l Res. 2005, 164, 811–828. [Google Scholar] [CrossRef]
  21. Kapitsa, S.P. Obshchaya teoriya rosta naseleniya Zemli. In The General Theory of World Population Growth; Nauka: Moscow, Russia, 1999. [Google Scholar]
  22. Popkov, Y.S. Mathematical Demoeconomy. Integrating Demographic and Economic Approach; De Gruyter: Berlin, Germany, 2014; p. 534.
  23. Kaashoek, M.A.; Seatzu, S.; van der Mee, C. Recent Advances in Operator Theory and Its Applications; Springer: Berlin, Germany; Heidelberg, Germany, 2006; p. 478. [Google Scholar]
  24. Rubinstein, R.Y.; Kroese, D.P. Simulation and the Monte Carlo Method; John Wiley & Sons: Chichester, UK, 2008. [Google Scholar]
  25. Popkov, Y.S. Macrosystems Theory and its Applications; Springer: New York, NY, USA, 1995; p. 340. [Google Scholar]
  26. Jaynes, E.T. Information Theory and Statistical Mechanics. Phys. Rev. 1957, 106, 620–630. [Google Scholar] [CrossRef]
  27. Popkov, Y.S.; Popkov, A.Y. New Methods of Entropy-Robust Estimation for Randomized Models under Limited Data. Entropy 2014, 16, 675–698. [Google Scholar] [CrossRef]
  28. Ioffe, A.D.; Tihomirov, V.M. Theory of Extremal Problems (Studies in Mathematics and Its Applications); Nauka: Moscow, Russia, 1979. [Google Scholar]
  29. Alekseev, V.M.; Tihomirov, V.M.; Fomin, S.V. Optimal Control; Nauka: Moscow, Russia, 1987. [Google Scholar]
  30. Darkhovskii, B.S.; Popkov, Y.S.; Popkov, A.Y. Monte Carlo Method of Batch Iterations: Probabilistic Characteristics. Autom. Remote Control 2015, 76, 775–784. [Google Scholar] [CrossRef]
  31. Gonzalo, J.A.; Munoz, F.-F.; Santos, D.J. Using a Rate Equation Approach to Model World Population Trends. Simul. Trans. Soc. Model. Simul. Int. 2013, 89, 192–198. [Google Scholar] [CrossRef]
  32. Shryock, H.S.; Siegel, J.S. The Methods and Materials of Demography; United Nations Department of Commerce: Washington, DC, USA, 1973. [Google Scholar]
  33. World Population Prospects: The 2015 Revision, Methodology of the United Nations Population Estimates and Projections; Working Paper No. ESA/P/WP.242; Department of Economic and Social Affairs, Population Division: New York, NY, USA, 2015.
  34. Shampine, L.F. MatLab Program for Quadrature in 2D. Appl. Math. Comput. 2008, 202, 266–274. [Google Scholar] [CrossRef]
  35. Coleman, T.F.; Li, Y. An Interior, Trust Region Approach for Nonlinear Minimization Subject to Bounds. SIAM J. Optim. 1996, 6, 418–445. [Google Scholar] [CrossRef]
Figure 1. The joint PDF of the random parameters r and u r for the range J = I r I u r .
Figure 1. The joint PDF of the random parameters r and u r for the range J = I r I u r .
Mathematics 04 00016 g001
Figure 2. The family of the PDFs of the noise ξ i , i [ 0 , 7 ] .
Figure 2. The family of the PDFs of the noise ξ i , i [ 0 , 7 ] .
Mathematics 04 00016 g002
Figure 3. The ensemble of projection trajectories on the testing interval.
Figure 3. The ensemble of projection trajectories on the testing interval.
Mathematics 04 00016 g003
Figure 4. Box plot for the ensemble of projection trajectories on the testing interval.
Figure 4. Box plot for the ensemble of projection trajectories on the testing interval.
Mathematics 04 00016 g004
Figure 5. The ensemble of projection trajectories on the prediction interval.
Figure 5. The ensemble of projection trajectories on the prediction interval.
Mathematics 04 00016 g005
Figure 6. Box plot for the ensemble of projection trajectories on the testing interval.
Figure 6. Box plot for the ensemble of projection trajectories on the testing interval.
Mathematics 04 00016 g006
Table 1. Estimation interval T e s t .
Table 1. Estimation interval T e s t .
i01234567
year19601965197019751980198519901995
E real est 3.0263.3583.6914.0704.4494.8845.3205.724
Table 2. Testing interval T t s t .
Table 2. Testing interval T t s t .
i01234
year19952000200520102015
E real tst 5.7246.1286.5146.9167.359
E 1985 prn 5.6665.9626.4506.9857.469
Table 3. Forecasting interval T f r c
Table 3. Forecasting interval T f r c
i012345
year201520202025203020402050
E UN frc 7.3597.6447.9648.2848.9249.564
Table 4. Calculated Lagrange multipliers.
Table 4. Calculated Lagrange multipliers.
Range I r , I u r
Measurements01234567
θ ¯ 0.0000−0.3833−0.3984−0.5839−0.3802−0.4679−0.18120.8881

Share and Cite

MDPI and ACS Style

Popkov, Y.S.; Dubnov, Y.A.; Popkov, A.Y. New Method of Randomized Forecasting Using Entropy-Robust Estimation: Application to the World Population Prediction. Mathematics 2016, 4, 16. https://doi.org/10.3390/math4010016

AMA Style

Popkov YS, Dubnov YA, Popkov AY. New Method of Randomized Forecasting Using Entropy-Robust Estimation: Application to the World Population Prediction. Mathematics. 2016; 4(1):16. https://doi.org/10.3390/math4010016

Chicago/Turabian Style

Popkov, Yuri S., Yuri A. Dubnov, and Alexey Yu. Popkov. 2016. "New Method of Randomized Forecasting Using Entropy-Robust Estimation: Application to the World Population Prediction" Mathematics 4, no. 1: 16. https://doi.org/10.3390/math4010016

APA Style

Popkov, Y. S., Dubnov, Y. A., & Popkov, A. Y. (2016). New Method of Randomized Forecasting Using Entropy-Robust Estimation: Application to the World Population Prediction. Mathematics, 4(1), 16. https://doi.org/10.3390/math4010016

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop