- Email: [email protected]
The parameter identification of a population model of China
0005 - 1098/84/$ 3.00 + 0.00 Pergamon Press Ltd. ~ 1984 International Federation of Automatic Control
Automatica. Vol. 20, NO. 4. pp, 415-428. 1984 Printed in Great Brilain.
The Parameter Identification of a Population Model of China* YUE-XIN Z H U t and BAI-WU WAN++
A mathematical form for the survival function drastically reduces the degrees of freedom in the existing population model and makes the solution of this parameter identification problem, under very limited statistical data, practicable. Key Words--ldentification; least-squares estimation; modelling; optimization; parameter estimation; population model.
death population in several successive years to estimate the age-specific mortalities directly (Boyarsky and Shusherin, 1955; Keyfitz, 1968). But these data do not exist in the ordinary census in China. That is why population model parameters for China have not yet been obtained. Hence, the study of estimating the population model parameters on the basis of the statistical data currently available is important. Our goal is, on the basis of the statistical data available, which include the data from censuses in the past 25 years, to study the parameter estimation problem for a population model of China by system identification, and to try to acciuire some useful model parameters. By comparison with common physical systems, one of the features of population model parameter estimation is that past statistical data, no matter how limited, are the only window from which the behaviour of the population process can be observed. The statistical data from the Chinese population system which are available for parameter estimation are very limited. Another feature is that the data-sampling time-periods of different output data of the population system are not identical, so the model structure assumes a very complicated form. Thirdly, the population model is a bilinear control system model with a rather high dimensionality and quite a lot of parameters to be estimated. Finally, the irreversibility of the population process means that the population forecast and, therefore, the prediction of model parameters themselves, have great importance. All these features give rise to the difficulties in the parameter estimation of a population model of China. In this paper five steps are given: (1) we choose a discrete mathematical model to describe the population process; (2) outline the available statistical data and their preprocessing; (3) transform the common model structure into a form
A l ~ ' a e t m A mathematical model for predicting and controlling the future trends of the population of China is given. The model used is a discrete bilinear system with 101 dimensions. The model is decomposed, transformed and its parameters are estimated by nonlinear optimization and curve-fitting techniques. Results for both time-invariant and time-varying parameter estimation are obtiined. Finally, a successful validation of the population prediction over the past 25 years is given, in which the predicting accuracy ofthe total population for each year is within about 1 ~. This is despite the fact that the available statistical data for these years, which include two censuses, are limited, incomplete and contain numerous contradictions. INTRODUCTION
CONTROL theory has been finding important applications in the study of the population process. With the deepening of these quantitative studies, one problem in particular has arisen--the parameter estimation problem of a population model, because it directly affects the credibility of the process forecast, the reliability of its synthesis and, therefore, population policy decisions. However, the population model parameters, i.e. the age-specific survivals and the age-specific fertility rates, have all traditionally been evaluated by means of demographic methods which require much more detailed statistical data than ordinary census data. For example, those demographic methods need the statistical data for the age-specific distribution of
* Received 30 March 1982; revised 2 December 1982; revised 13 October 1983. The original version of this paper was presented at the 6th IFAC Symposium on Identification and System Parameter Estimation which was held in Washington, D.C. during June 1982. The published proceedings of this IFAC meeting may be ordered from Pergamon Press Ltd, Hcadington Hill Hall, Oxford OX3 0BW, U.K. This paper was recommended for publication in revised form by associate editor J. Mendel under the direction of editor H. Kwakernaak. t Department of Electric Engineering, Xi'an Institute of Optics and Precision Mechanics of Academia Sinica, Xi'an, Shaanxi, People's Republic of China. :~Institute of Systems Engineering, Xi'an Jiaotong University, Xi'an, Shaanxi, People's Republic of China. 415
416
YUE-×IN ZHU and BAI-WU WAN
corresponding to the real data structure; (4) decompose it further into two subsystems--the age-process and t he birt h-process, and identify them separately to simplify the problem; and (5) formulate the parameter estimation problem and convert it into an optimization problem with dimension 101. In spite of this, it is still too difficult to estimate the model parameters of either the age-process or the birth-process directly because the statistical data are too limited. So we find some relations among the components of the parameter vectors under estimation by a curve-fitting technique, and thus simplify the nonlinear optimization problem by reducing its dimensionality. We finally show that the optimal time-varying parameter estimates, as well as the time-invariant ones, can both be found, and, give a satisfactory validation. MODEL STRUCTURE AND ITS TRANSFORMATION
Choice o f the mathematical model
Several kinds of mathematical model have been developed to describe the population process. They can be briefly summed up as follows: first, continuous mathematical models with continuoustime, continuous-age and continuously-distributed state variable, which are expressed as partial differential equations (Foerster, 1959; Song and Yu, 198 l). Second, discrete models with discrete-time and discrete-age variables and lumped-state variables, which are mainly expressed as a Leslie matrix or one of its variants (Leslie, 1945; Kwakernaak, 1977; Wang, 1980). Third, a model characterized by a finite-state continuous-time semi-Markov process (AIj and Haurie, 1980). Given a large scale of China's population system and its statistical data structure, we prefer a discrete model to describe the population process. The general structure of this model is (Kwakernaak. 1977; Wang, 1980) x(k + 1) = H ( k ) x ( k ) + qo~(k)B(k)x(k) + w(k) (1) where k is a discrete-time variable, in year, with origin at 0AD (the moment k, refers to the appointed data-statistical time-moment in the chosen year kAD)
a discrete-age variable. The oldest age (considered in the model) is m = 100
H(k) =
.
i]
~7~(k) ....
(2)
r/,._ l(k)
is a Leslie matrix with dimension (m + l )2, in which all the elements are equal to zero except t/a(k), a = 0, 1 ..... m - 1, on the first lower-left subdiagonal of it. In (2) G ( k ) = xa ÷ i (k + 1 )/x,(k),
a = O, 1. . . . . m - 1
(3a) ~loa(k) = xo(k + l)/'yo(k) = r/_l(k)
(3b)
are called the survival at age a and the baby survival in the very birth-year (Wang, 1980), respectively, at the moment k, in which yo(k) is the aggregate number of babies born in [k, k + 1). For simplicity, we define t/'(k) = [t/oa(k) r/o(k)..... r/,._l(k)]
(4)
as the age-specific survival vector with dimension m + 1. And B(k) = [b(k) i O(m+l)×m]'
is the birth-matrix where b'(k) = [bo(k)
b l ( k ) ..... b,.(k)]
(5)
is the age-specific fertility-rate vector and ~0, O > a o r a > L bo(k) = [yo.Ak~/x~(k), 0 < a < ~
(6)
where yo.~(k) represents the number of babies born by x~(k) in [k,k + I); and 0 and ~ are the youngest and the oldest fertile-ages, respectively. In this paper, 0 = 17, ~ = 50. The vectors q(k) and b(k) are the model parameters to be estimated. The term b(k) is also usually written as b(k) = k:(k)fl(k);.(k)
(7)
where k : ( k ), fi{k) and ).(k) are the ratio-coefficient of women in fertile-population, the average fertilityrates of women and the fertility-pattern vector, respectively"; and
x'(k) = [xo(k)xt(k) . . . . . .'¢~(k)] (8a)
fl(k) = V h,,(k),k/.(k}
and
a =0
2'(k) = [2o(k)2~(k) ..... 2~(k)l
w'(k) = [wolk)wt(k) . . . . . w.,(k)] are, respectively, the population state vector and the population state disturbance vector at the moment k, the prime indicates a vector transposition; x~lk) and wAk), a - 0, l ..... m, are the aggregate number of people and the number of disturbance people, respectively, in age-region ,~a.a + I ). The a denotes
2,,(k) = b,Ak)/{kr(k)fl(k)~.,
(8b)
a = O. 1. . . . . m
(8c) and m
Y :.,,(k)= 1.
The parameter identification of a population model of China Obviously, the birth-process can be more clearly controlled and described, and therefore identified in terms of it(k) and fl(k) than in terms of b(k).
417
z',r(k) = [yo(k)yo(k + 1). . . . . yo(k + n ,
1) 3 (16)
and [
z~o(k + n) = [x'(k+ n) , s;(k) ',,d'.(k)] Transformation of the model structure Equation (1) is not suitable for the parameter estimation because it is difficult to identify this bilinear system directly; its discrete-time period is different from the real data-sampling time-periods; and there is no output equation corresponding to real statistical data. Hence, the model structure (1) must be transformed. We first divide the model (1) into two subsystems: the age-process subsystem x(k + 1) = H(k)x(k) + r/oQ(k)yo(k) + w(k) (9) and the birth-process subsystem yo(k) = B(k)x(k)
(10)
y~(k) = [yo(k) 0,..., 0 ]
(11)
where
is the birth-vector with dimension 1 x (m + 1). These two subsystems are common linear systems, directly coupled with each other, and both are concerned with only one of the two parameter vectors ~/(k)and b(k), respectively. So we can identify them both separately to avoid identifying the bilinear system directly (Eykhoff, 1974; Isermann, 1979). Next, assuming the data-sampling time-period about the population states to be n years, we derive the model structure corresponding to the data structure, including the output equations, as follows:
are the output vectors of the birth-process and the age-process, with dimensions 1 x n and 1 x (m + 2n + 1), respectively, s',(k) - [s(k + 1)s(k + 2). . . . . s(k + n)]
[d(k)d(k + l) ..... d(k + n - 1)] (19)
d,,(k)
are the total population vector and the total death population vector in [k, k + n), respectively, where s(k + j) and d(k + j - 1),j = 1, 2 . . . . . n, are the total population at k + j and the total death population in [k + j - 1, k + j), respectively. The parameter matrices in (12) and (13), can be briefly written as h
Hha,(k) = l-I H(k + p - j ) ,
(20a) in particular
Hp,r(k ) = Hp(k),
H, ( k ) = H(k),
where I,"+, is the identity matrix with dimension (m + 1)2; and
Hi(k) = dia
--ri
L.i=
th(k + h + i)
- 1
0 0
17[ th(k + h + i - 1)..... i= -1
1
(12b)
for the age-process subsystem, and yo,(k) = B,(k)x(k) + B~,(k)w,,(k)
(13a)
z~r(k) = Eyo,(k)
(13b)
h_2 I-I rh(k + i +
1
1) 0 , . . . , 0 , h = I, 2 . . . . ,n
i=-1
h-
1
h,...,m
for the birth-process subsystem where
with dimension (m + 1)2; and
y~,(k) = [yo(k + n - 1)yo(k + n - 2).....
H~(k) = [Ho.h(k) ,[ H, .h(k)', .....
(21)
' [Hha.h(k)], h = 1, 2 . . . . . n
(14)
(22a)
and
and w'.(k) = [w'(k + n - 1) ,l w ' ( k + n - 2 ) ]
(20b)
Ho.p(k) = Ho(k)= 1,,+1
z,,,,(k + n) = Co(k)x(k) + Do(k)yo,,(k)
yo(k)O ..... 0]
h = 1,2,...,p ~< n
j=l
(12a)
+ Ro(k)w,,(k)
(18)
and
x(k + n) = H,,(k)x(k) + H~(k)yo.(k)
+ H~,(k)w,,(k)
(17)
! h H~l(k) = [0(,"÷1)x.,-h),Hw(k)]. h = 1,2 . . . . . n
, ..... w'(k)]
(15)
are the total birth-vector and the total disturbance vector in [k, k + n), with dimensions 1 x (m + l) and 1 x n(m + 1). respectively
(22b) with dimensions (m + l) x h(m + l) and (m + 1) x n(m + 1). respectively, and, when h = n
H~, (k) = H~(k)
(22c)
418
YUE-X~N ZHU and BAz-WU WAN I
B'.,(k) = [H~,_t(k)b(k + n - l) I
/{oo(k) = diag[l - qoo(k + n - 1) t
t
1 - qoo(k + n - 2) .....
H/,_2(k)b(k + n - 2)~r.... i I !
I ))"
H;(k)b(k))O~
..,i 0]
1 - qoo(k) 0 . . . . . 0]
(23)
and
(27c)
both with dimension (m + I)2; and
B~,(k) = [ H ~ ( ~ ( k ) b ( k + n - 1) . . . . . i
H~lt(k)b(k + 1)
0
r=
i . . . . . 'n, O ]
(24)
with dimensions (m + 1)2 and n(m + 1) x (m + 1), respectively. The observing matrices are (25a)
C'.(k) = [H,,(k) l f .(k) I FS(k)
[H,(k)~G~(k)~Gi(k)]
D'o(k) =
(25b)
both with dimension (m + 1) x (m + 2n + 1), and R'o(k) = [I-t"~ (k) l L,(k) l L'd(k) ]
(25c)
and
L1' ,1',, [o 0,0,×,,,_.,-1) 1
E =
(25d)
with dimensions n(m + 1) x (m + 2n + 1) and n x (m + 1), respectively. In the above, iI , l F'~(k) = [ n ' , ( k ) l Hz(k)l ~). . . . ..' ttH;(k)l ] d
(26a) i
F'n(k) = [H'o(k)l~(k) H~(k)l~(k + l)! ..... i H~_l(k)l~(k + n - 1)]
(26b)
G'~(k) = [H,~l(k)l ,')H2t(k)l iI , .... iIH,. l(k)l ]
(26c)
G'.,(k) = [0 H ~ ( k ) l ~ ( k + 1), .... , H~L-l(k)l~(k + n - 1)] + I~o~(k)E'
(26d)
all with dimension (m + 1) x n, and I
,2
I
I
t'~(k) = [H~ll(k)l ~H~t(k)l I ..... ~H~,(k)l]
(26e)
..,i i
I ~L L'~(k) = [O[Hwl(k)~(k + 1 ) II, .
,
' " - ' (k)l~(k - + n Hw~
1)]
(26f)
both with dimension n(m + 1) x n. Also in the above #'(k + j ) = [1 - qo(k + j ) 1 - r/t(k + j ) .....
j = 0 , 1 ..... n - 1 (27a)
l-qm_~(k+j)l],
with dimension 1 x (m + 1); and
[1
1,...,1]
the dimension of which corresponds to that of a matrix multiplied by the vector, 1, or that of another vector multiplying the vector 1 in equations. Equations (12) and (13) are typical and canonical linear control systems but their internal structure is rather complicated. Their derivation is given briefly in Appendix A. Matrices H,(k), H~(k), H~(k), Co(k), Do(k) and Ra(k) in (12) contain only t/, while matrices B,(k) and B~,(k) in (13) contain both q and b. So the age-process (12) can be identified solely by the parameter vector q alone. Hence, we can identify the birth-process after having identified the ageprocess. TRANSFORMING PARAMETER ESTIMATION PROBLEMS INTO OPTIMIZATION PROBLEMS Available statistical data and their preprocessing The original statistical data available for identification are not only limited and incomplete but also contain numerous contradictions. Therefore, we had to preprocess them by: standardizing, their data-collecting time-moment; converting the regional data into the national ones; converting and fitting the baby mortality data; investigating the precision of all data used; and further processing the parts of those data containing obvious contradictions and serious errors, and so on. Then we obtained all the preprocessed data, and their precision distribution, which consist of: (a) ilk), k = 1953, 1964, 1975, where the symbol " denotes statistics of the vector under it, of which i(1953) and i(1964) come from the 1953 and 1964 national census data, respectively; while i(1975) comes from the regional census data of 1975 and have been converted into national ones. Thus, there are only three sets of sampling data of the population states at three different moments in two data-sampling time-periods, both with n = 11 years. (b) g = [stk)s(k + 1)..... s(k + 25)] a ' = f l l k ) d ( k + 1)..... d(k + 24)] ~' = [Yolk)Yolk + 1). . . . . yolk + 24)]
(28)
-1
~~0, 1-I r/ilk + h + i) l ..... i
H i t ( k ) = diag E0
t=--
0..... n-h-1 n-2
i
n -
t
y~ ~,(k + i + 1)~0~ . . . . . " ',"
i=-1 n -
and
I
n
1
,...,
I
m
h
...
where k = 1953. These all come from routine annual statistical data, and their data-sampling time-period is one year. The formulation o f the parameter estimation problem
h = l ..... n (27b)
From (17), we see that the state vector x(k + n) is contained within the output vector z.o(k + n) of the age-process (12), which has been arranged on
The parameter identification of a population model of China purpose, and therefore, the state equation (12a) is contained within its output equation (12b), too. Hence it ~s enough for the age-process (12) to be studied through only its output equation (12b). Considering the recurrence properties of discrete models, we can represent the age-process as the model with input vectors x(k), yo.(k) and w.(k), and the output vector z..(k + n), shown in the dotdashed block in Fig. l(a), if there were no symbols ^ and - above the vectors. As for the birth-process (13), it is really a set of algebraic equations. Because of the special structure of matrix E, there is essentially no difference between the two vector, yo.(k) and z.v(k ). Hence, the birth-process can also be represented as a model with input vectors, x(k) and w.(k), and output vector yo.(k), which is shown in the dot-dashed block in Fig. l(b). Now we compare the mathematical model with its own real process for the age-process and the birth-process, respectively. As shown in Fig. l(a) and (b), respectively, we take the input vectors of a real process corrupted by measurement noises (i.e. statistical errors) as the input vectors of its model, and then use the output response of the model and the output vector of the real process corrupted by measurement noises to form the error vector as
er(k) = ~o,,(k) - ~o,(k) = B,(k)~(k) + B~,(k)f%(k) - ~o,(k)
Transforming the parameter estimation problems into optimization problems In the case of the age-process, considering the age-specific survival vector # as time-invariant during every time-period, and in the sense of weighted least squares, we transform the parameter estimation problem of the age-process into a nonlinear optimization problem, for every timeperiod respectively, as follows:
Co(k)f~(k) + D,,(k)~o,,(k)
(29)
+ Ro(k)~k,,(k) - ~,,°(k + n)
with dimension ( m + 2 n + l ) x 1 for the ageprocess, as shown in Fig. 1(a), or as
vz(k+n) w n(k) x(k)
[
Yon(k)
/
,.(,)
I
Zno( k +n J
,+++Z~lk+n)
't '1"" Ji;rm(k + n)
-leo(k+n)
ReoL age -process
,
I
Colk)
Do(k) Math. model of the age prooess
wn(k)
X(~)
,
'* '1 Ro( , •
]
[
ReaL
(a') _J
~.(k}
,.,,,
II I
L
i
vy(k)
~(~+ ;o. (k)
.
birth - process
;r =[ .:(k, I+.°o,+ v,lk),~÷ x(k)
Yon(k)
~l-*y (k)
(b) Math. rn~et of birth-process ----.-J
FIG. 1. The formulation of identification problems: (a) for the agc-proo:ss, (b)for the birth-process.
AUTO20:6-'C
(30)
with dimention (m + l) x 1 for the birth-process, as shown in Fig. l(b). In the above two equations and in Fig. 1, the overbar denotes a response of a model to its inputs. The vectors v=(k), vy(k), vw(k) and v.(k + n) are the corresponding measurement noises. Note that the error vectors, (29) and (30), are both expressed as functions of the parameter vectors r/ and b under estimation, respectively. What either of the two parameter vectors is taken as determines the corresponding error vector and, so determines the closeness between the model and its real process and, therefore, the effectiveness of the model, and vice versa. Hence, the specific formulation of this parameter estimation problem here ought to be: choose the model parameter vector so that the corresponding error vector, in a sense, is a minimum. It is worth noticing that these two parameter estimation problems are both problems in which the input a.nd output-signals are corrupted by measurement noises.
eo(k + n) = ~,(k + n) - $,,(k + n) =
419
420
YUE-XIN ZHu and BAI-WU WAN b'(k + j ) = [bo(k +j)bo+,(k +j) ..... bdk + j ) ]
,.a~-,
+ 2n + 1 e.(k +
(31)
where R" + 1 is a real vector space with dimension m + 1; J~(q) the objective function with the same dimension as the mean square error; e,(k + n, q)just the error vector e.(k + n) which is a function of q; Q the weighting diagonal matrix determined by the precision distribution of corresponding data. Thus the parameter estimation problem of the age-process has been transformed into a feasible but rather complicated and intrinsically nonlinear optimization problem with dimension 101. As for the birth-process, since b(k) appears linearly in B,,(k) and B~,(k), and also in the error vector ey(k), as seen from (23), (24) and (30), where q(k) has been considered as known in identifying the birth-process as mentioned above. So the parameter estimation problem of the birth-process is really the unique regression-solution problem of the following linear matrix equation (Eykhoff, 1974)
$o.(k) = B.(k, b),t(k) + B~(k, b)#.(k)
(32)
with respect to b, where B,.(k, b) and B~(k, b) are just B,,(k) and B~(k), respectively, and are linear functions of b(k). Since q is known, (32) can be written as
~o,(k) = [x'(k + n - l)b(k + n - 1},...I
01 where x(k + j ) , j = 1,2 ..... n - 1, can be derived from the real statistical vectors :~(k), w.(k) and the known vector q according to (12a). Considering the limitation of the birth-statistical data, we identify the birth-process not only within a s t a t e - - d a t a sampling time-period but all over the 25 years for which we have birth-statistical data, as shown in (28). And, for simplicity, when treating the birthprocess, we cut out the redundant zero elements in vectors yo.(k), b(k) and A(k), and represent them with the same symbols as above. Thus, the parameter estimation problem of the birth-process can be finally transformed into the unique regression-solution problem of the linear matrix equation, either ~ ' = [x;(k + l -
l)b(k + l -
and x;,(k + j ) = [xo(k +j) Xo+,(k +j) .....
x;(k+j)],
j = 0 , 1 ..... l -
both with dimension 1 x ( ( - 0 matrix X ' = [xb(k + l -
1
(35)
+ 1); and the
l)!xb(k + l - 2) i... I Xb(k)]) I (36)
with dimension (~ - 0 + 1) x I. In the same sense as the above, xb (k + j) and therefore X are considered as known here. IDENTIFICATION OF THE AGE-PROCESS
The solution of the optimization problem (31) and difficulty with the parameter estimation In order to solve the optimization problem (31), we used a TQ-6 computer (a large computer made in China); and used the Davidon-Fletcher-Powell optimizing algorithm which requires an analytical expression for the gradient of the objective function. We solved the optimization problem for the two time-periods, i.e. 1953-1964 and 1964-1975, respectively, and then obtained two optimal estimates of the age-specific survival vector. A typical example of these two results, for 1953-1964, is shown by the discrete-point set in Fig. 2. The other results are similar. Obviously, such a result is impractical. In other words, these optimal parameter estimates, from a direct solution of (31), are too dispersive either to be taken as the final result of the parameter estimation, or to be further processed by means of the usual techniques of graduation of data. This is just the difficulty with this parameter estimation by a direct solution of optimization problem (31). Obviously, in view of system identification, that is because of too limited statistical data.
Analytic relation among the age.specific survivals Concept. In the above as well as usual discrete population models, all the components of the vector ~/have been regarded as independent. However, they are interrelated through the human-fife-process in a
l)i.., i i¢;(k)b(k)] i
°
oI
(33)
:
with respect to b(k + j),j = 0, l . . . . . l - 1, which are assumed to be time-varying; or 9
= Xb
(34)
"~. oaok I
with respect to b which is assumed to be timeinvariant. In (33) and (34), k = 1953; l = 25 denotes the number of the birth-statistical data used ; ~ is just the birth-statistical data vector (28), with dimension lxl
" ""
070~-
.
• \
i
~°6° I-F-
:~_ 050
0
.,~
~
r
I
I
I
1
I
I
t
tO
20
30
40
50
60
?0
BO
90
Age
(yeors)
FIG. 2. A direct solution of (31).
,
The parameter identification of a population model of China society so that among them there would exist a relation which is the very thing we are interested in. On the other hand, strictly speaking, the survival in this paper as well as in demography is, in nature, a probability of survival. However, the surxlval defined in (3), i.e. the survival quotient for individuals, is only a frequency of occurrence of survival under the condition of finite individuals. Obviously, the probability of survival is determined by social factors, medical conditions and so on, while the closeness of the frequency to its probability is determined by the number of individuals under discussion. The relation among the age-specific survivals discussed here is just the dependent relation among the individual probabilities of survival for the different age-set population. Hence, this relation is deterministic. For simplicity, we will still use (3) as the definition of the age-specific survivals, but understand it in the same way as the above.
The survival function. On the basis of continuous population models existing, we now have a proposition as follows. Proposition. For a set of the age-specific survivals in a discrete model at a moment k, r/o(k), a = - 1, 0 ..... m - 1 , there must exist a continuous and differentiable function r/(~), with respect to the continuous age-variable ~, which, in the sense ofthis discrete model, is just equal to the age-specific survival at the age ~ e [ - 1, m - 1 ] and at the same moment k. In particular, when ¢ = a, the function r/(¢) satisfies r/(a) = ~/,(k), a = - 1 , 0 . . . . . m - 1.
421
and invasion of diseases, reach some equilibrium, and the survival reaches a maximum. After that, the increase in the factors resulting in death exceeds the influence of the increase in adaptability, and the survival decreases slowly over a rather wide agerange. After about 50 years of age, for example, senescence begins to increase, and after some age, for example 60, the survival starts decreasing quickly. The greater the age, the faster the decrease of the survival until the survival vanishes. That is the rough shape of the survival function, and therefore, its degree of freedom has been restrained.
An analytic relation. From the above, the teenager section of the survival function can be approximated by an exponential function; the middle-age section by a linear function; and the old-age section by a composite power and exponential function. Summing these, we then obtained an approximate expression for the survival function with seven independent parameters
-
-
ctaexp(--at2~ ) +
expt - 0,• I ~e [-1, m)
) (38)
where
is the independent parameter vector with dimension 7; the factor in the linear term
(37)
The proof of the proposition is given in Appendix B. Obviously, the function ~/(¢) is just the relation we are looking for, and we call it the age-specific survival function, or simply, the survival function. In order to find the specific expression of this function, let us look at the real statistical law of the age-specific survivals.
Statistical laws. Many real statistical results show that the age-distribution of the survival has the following laws: The survival of a group of new-born babies is generally lower because their existing state changes rapidly and they are poorly adapted to their new surroundings. As their age increases, their adaptability increases quickly and so does their survival. However, the older they are, the lower the rate of increase of their survival becomes. On the other hand, with the increasing age, the probability that diseases and other external harmful factors will affect people rises slowly but surely, and almost linearly. At some age in the teens, the two factors, the increase of adaptability to external surroundings
J is chosen so that it does not greatly influence the value of ~(~, e) when ~ < 70, but it gives ~(~, ~) the following analytic features
> 0,
[ - 1 , oo)
and as ~ --* oo, r/(~, at) --. 0 monotonically.
Effect. Is it effective that the survival function (38) is used to express the age-distribution of the survival? Obviously, this can not be proved theoretically and can only be verified by practical survival data. An effective survival function should be able to approximate various kinds of age-specific survivals closely enough, by means of changing the independent parameters in it. In order to verify the function (38), we used it to fit several kinds of agespecific survival data which came from real, published statistical results and had been graduated. This is a curve-fitting problem, i.e. the solution of the nonlinear optimization problem
422
YUE-XIN ZHU and BAI-WUWAN min J,(,,) =
where R 7
/ [ m - - ~ ],q(at) - @H~1
min Jo{q(=))
(39)
~t~R 7
,
is real vector space with dimension 7
~'(,')= [ q ( - l , ~ )
q(O,~)...q(m- 1,=)]
(40)
is the fitting function vector; and
#'= [qo,, @o.-.0,.-i]
(41)
is the age-specific survival data to be fitted as a vector. This was done with the DFP-algorithm on a DJS-130 computer (a small computer made in China). The results of fitting this function to many sets of real statistical data show that the survival function (38) is satisfactory. As an example, Fig. 3 shows a result of fitting (38) to a kind of the age-specific survival data, which are published in (Song and Yu, 1981). The largest fitting errors are within I%0, when age <57; and within 1 ~ , when age <70. And the fitting errors to the individual probabilities of those age-specific survival data must be much less than the above. Hence, function (38) is the relation among the components of the age-specific survival vector we are looking for.
Practical identification algorithms Using the curve-fitting method of (39), agespecific survival data like the direct solution of (31) can be smoothed easily no matter how dispersive they are, and therefore the parameter estimation problem of the age-process can be solved. This, however, requires two steps: first, solve (31), and second, solve (39). Now substituting p/(at) in (40) for q in (1), we obtain a new model structure of the age-process which has only seven independent parameters. Substituting the q(=) for t/in (31), we can therefore rewrite the parameter estimation problem of the age-process as LefL scare .u U
Time-invariant and time-varying parameter estimation results of the age-process Assuming the age-specific survivals to be timeinvariant in the two periods, 1953-1964 and 1964-1975, respectively, and following the above method, we then obtain two sets of optimal parameter estimates of the age-process, q,(~) and r/z(~), which we call the time-invariant results of identification of the age-process. In order to obtain the dominant law of the agespecific survivals with respect to time, we noted the following facts: first, the total mortality of the Chinese population has, basically, been decreasing exponentially over the 25 years, with the exception of a few years of unusual mortality--the unusual part of which has been considered as an external disturbance; second, in a steadily developing society and under the condition that the senectitude states are basically stationary during the time-region under discussion a reduction in the total mortality must, to some extent, correspond to some reduction Rig.t
I
,...
r
scale
0 80
0 99
i, I,
I
"%~
Fitting 098 • o,
Fi%tea
t
060 ~-
o,~ e,,,. 0
2,0
40
60
(42)
instead of (31). It can be seen that the dimension of the nonlinear optimization problem is reduced from I01 to 7, and its computation-time is evidently reduced too. Using the DFP-algorithm without the analytic gradient, we can also do it on a DJS-130 computer. Since relation (38) represents the internal relation among the components of the age-specific survival vector q, it can be regarded as a further identification of the usual discrete model structure, and this is just the distinctive point of this paper in treating the identification problem. And that is why we call this kind of parameter estimation a parameter identification.
I
m
•
q
+ 2n + t Ile,{k + n, 1/(~t)}lt
80 Age (years)
I00
Fla. 3. A real example of fitting (381.
{20
140
The parameter identification of a population model of China of the age-specific mortalities. So we can assume the dominant law of the age-specific mortalities with respect to time to be exponentially decreasing, too; and ~h (~) and r/2(¢) obtained above to be equal to the age-specific survivals at the moments t~ and t2 = tl + n, respectively. Thus the general expression of the survival function can be written as r/(¢, t) = 'l ~.(¢) - c(~) exp { - :t(¢)t } ,1(¢.',tl ) = '11 (¢') ,7(¢,
(43)
t2) = ~2(¢)
where q=(~), c(~) and ~(~) are functions to be defined. Furthermore, we assume ~/=(~.) to be dependent on q:(~) and find the following form of ~(¢, t)
,ff~, t)
=
1 - k; + kcr/:(~) - [1 - k; + keq2(¢) - )h(~.)] 1 (1 -
k:)[1
-
"~.,+,)/.
r/=(¢)]
-
(44)
where k{ and t~ are both parameters to be determined under a definite time-origin. Thus the time-varying parameter estimation problem of the age-process can be roughly formulated as: choose the parameters k~ and t] so that the precision of total population forecast for each year in the past 25 years is optimal according to the survival parameters determined by (44). Thus we obtain the time-varying results of identification of the age-process which are expressed by the survival function (44) with respect to ~ and t, by means of which we can evaluate the age-specific survivals at any moment, not only from 1953 to 1978 but also after 1978. The former are called timevarying identification results and the latter the I /
parameter forecasting (or extrapolating) results corresponding to the former. When new statistical data come, we can use them to correct the earlier results, to obtain new results, and to make them more precise.
Validation of the identification results of the ageprocess Validity of the identification results can be verified only by studying the forecast precision of the past behaviours of the process, i.e. taking the population states in 1953 and 1964, respectively, as known initial distributions; the birth population every year as the inputs; and the parameters obtained above, time-invariant and time-varying, as those of the process model. We have calculated the population states in 1975, the total population and the total death population for each year in the 25 years, and compared them with real statistical data in Figs 4-6. Validation of the results with respect to both the time-invariant and the time-varying parameters obtained is basically satisfactory. The forecast precision of the total population for each year over the 25 years is within about 1 ~ . The validity and the extrapolation properties of the time-varying results are better than those of the time-invariant ones. So the forecast results of the parameters are valid too. IDENTIFICATION OF THE BIRTH-PROCESS
As mentioned in Appendix C, the identification of the birth-process here depends on the determinability of parameter estimates, or rather, it depends on whether the regression-solution of (33) or (34) is unique or not. Hence, if the time-varying parameter vector b(k) is considered as the direct object to be estimated, (33) is obviously not determinable. Even
J95~ -!975
:o
:
/~
v.:-
I
1964-1975 ~
°
k/
i
:o /),o ,i ,'o
'h//"l,"l
',1 'h
J964- 1975~
("~
(oI
1953-1975
~.ge (yeors ]
e,, -:0
423
/'~
I ' ~ I' " # I. I f % I/
(b)
FIG. 4. The predicting precision of the states in 1975: (a)by time-invariant results, (b)by time-varying results.
424
YUE-XIN ZHU and BAI-WU WaN
eee
*,eo
} Trine-varying result
...... I Time lrlvor~ont • • • • • A~ r e s u l t -
+I* Q
3
:
s
Q)
-S
t ,
\
io • ~ 3 ~ ~0
~58
,
• •~ • ~A73
,
__.&,
-,5 -- -
1953-1978
*
FIG. 5. The predicting precision of total population.
L 40
__-~____ } 1 9 5 3 - 1 9 7 8 • • • •
IA
r,'\--\
oI [ -'°b
',
',
..,.
..,-.m:" .o.
I ~.i¢ A" I *e~-'k I ~e 63T i ee r3o%~,,le78 ..... ........ " ' -
-e°/~-
. . . . ~ Time-varying • • • • • JI r e s u l t ...... Time-invoriant • • A• • result
F I G . 6. T h e p r e d i c t i n g
precision
of the total death
population.
the birth-process is mainly determined by the rank of X (see Appendix C). Obviously, if the population state is stationary (Kwakernaak, 1977), then, no matter how numerous the birth-statistical data are, the birthprocess is still not determinable. Although a theoretically stationary population does not really exist, a population process close to a stationary one, is generally not determinable either, because it usually requires finding the inverse of an illconditioned matrix. Thus, it follows that this kind of data-structure of the Chinese population system does not have universality in identifying the birthprocess, although it can be improved, to some extent, by means of a special algorithm. Fortunately, population states of China were not stationary in the 25 years and the rank of X, at least, is more than 10 which can be seen from its population states. In the Chinese population model, b(k) changes with respect to time mainly due to fl(k), while the fertility-pattern vector 2(k) is usually rather stable. So we use the following smooth function
b(~, t, v) = k/fll (t, v)21 (~, v) if b is assumed to be time-invariant, (34), and therefore the birth-process, is not determinable, either. This is because (rankX)m~=l<~-0+
1
where l = 25 as in (34); in other words, the birthstatistical data are too limited. However, as seen in identifying the age-process, the components of h are not independent of each other; a lot of real statistical results concerning the age-specific fertility-rates (for example, Ravenholt and Chao, 1974), show that there is some relation among them. In fact, nowadays some smooth functions are often used to express the dependent relation among the age-specific fertility-rates (Song and Yu, 1981). The most common one is the X2-probability density function
1 En/2- 1 e - ~ , 2 , 2,:2F(n/2) Px"(~) =
0
~< 0
e>0 (45)
where ~ = ~ - 0 + I; and n is a positive integer. This is used to describe the fertility-pattern function (Song and Yu, 1981), or rather, the relation among the components of the fertility-pattern vector 2, which is the same as b except for a scalar coefficient k:~ [see (8)]. On the other hand, it can be proved that if the relation among the components of b can be described by a smooth function with r independent parameters, and so long as rank X >/ r, (34) is then determinable. In the other words, with respect to our particular data structure, the determinability of
(46)
to identify the birth-process, in which k: is taken as a constant, equal to 0.4835 (Song and Yu, 1981); v' = [% vt...v9] is the independent parameter vector with dimension 10; while t - 8.005 fll(t,v) = Vo - vlt + v2tan - t - 0.001
exp{-O.5(t-v31V'~,t~[1953,1978], ~V,~l)
(47)
describes the time-varying properties of b and is determined by studying the change of the total birth-rate in the 25 years 21(¢, v) =
-
0 =, 1
'"
exp
-0.5
"V6
- 0 + 1
V8
(48) which describes the age-distribution of fertility-rate, or rather, the fertility-pattern, and is just a generalization of (45). In fact, it is the same as (45) provided that, in (48) v6=% =vo=
1 and
vT=n/'2- 1
except for the factor 1/[2"'2F(n/2)] which has been considered in (47). Thus, discretizing (48) and (47), according to (8) and (46) we can write the fertilitypattern vector 2' = [2o )-0~~... 2~] as
2 = 2t (v)/~/,.~ (a, v) where k't(v)= [).lt0, v)
,;.l(0+ 1, v)...21(~,v)]
(49)
The parameter identification of a population model of China
425
and the average fertility-rate of women as fl(k + j) = fl1(k + j,v) T. }q(a.v), a=O
j = O , 1..... 1 - 1.
(50) FIG. 7. The validation of the identification result of the birthpro~ss.
Substituting (49) and (50) into (33), we can write the linear regression equation (33) as the following nonlinear regression equation ~'
=
vector, obtained in this paper, as the model parameters, and by using the average fertility-rates of women, /3, as the parameter governed by the national family-planning policy. From Fig. 8, it follows that, in order to carry out the national population plan of China in which the total population is specified as under 1200 million by the end of this century, the average fertility-rates of women must be controlled under 1.70, and then the maximum of the total population will be about 1310 million and will take place in 2026 A.D.
DX).I (v)
where D = k/diag [fl~(k + 1 - 1, v) fl,(k + l -
2,v)... fll(k,v)].
So the parameter estimation problem of the birthprocess can be then converted into the nonlinear optimization problem ...,o
7 IIDX & (v) -
1
~'ll&.,
CONCLUSIONS
by which the birth-process parameter estimates, i.e. the fertility-pattern vector and the average fertilityrate of women, can be obtained. These results can be validated by the precision of calculations of the birth population in the 25 years, as shown in Fig. 7. Except for a few points, their precision is within 5 ~o-
APPLICATIONS
The parameter estimates obtained in this paper, along with the model, have been used to study prediction and control of the future trends of Chinese population in order to provide background material for the national population and familyplanning policy decisions. For example, Fig. 8 presents a few predictions for the total population, which were acquired by using the time-varying agespecific survival function and the fertility-pattern
From the above we conclude that the series of methods for identification seems to be effective in solving the parameter estimation problem of a population model under the conditions of very limited statistical data. This series includes implementation of a discrete model, decomposition and transformation of the bilinear model, a curve-fitting technique and nonlinear optimization, and, most important, the reduction of its dimensionality. The applicability condition for the curve-fitting technique is that there exists a functional relation among the parameters under estimation. The time-varying result makes the model more effective. This series of methods may be useful for population model identification in developing countries. Acknowledgements--The authors are indebted to Professor Bao~ sheng Hu, Huan-chen Wang, and Zhao-yong You for the discussions and their valuable suggestions.
t
B=2.20
g ,4 -._L3d_ _ _ . / ~ - e ~ - _ _ _
,o I.--
II
2026 6-
J980
2000
I
2020
I
2040
1
20¢~
I
~POQO
A.D. ( y e o r )
FIG. 8. A few predictions of total population of China by means of the identified parameters.
426
YUE-XIN ZHU and BAI-WU WAN
REFERENCES AIj, A. and A. Haurie (1980). Description and control of continuous-time population process. I E E E Trans. Aut. Control, AC-25, 3. Boyarsky, A. Ya and P. P. Shusherin (1955). Demographic Statistics. Gosstatizdat, Moscow (in Russian). Eykhoff, P. (1974). System Identification. John Wiley, London. Foerster, H. V. (1959). Some remarks on changing populations. In F. Stohman (Ed.), The Kinetics of Cellular Proliferation. Grune & Stratton, pp. 382-407. Isermann, R. (Ed.) (1979). Tutorials on system identification, 5th I F A C Symposium on Identification and System Parameter Estimation, Darmstadt. Keyfitz, N. (1968). Introduction to the Mathematics of Population. Addison-Wesley. Kwakernaak, H. (1977). Application of control theory to population policy. In A. Beusouson and J. Lions (Eds), New Trends in Systems Analysis. Springer, pp. 359-378. Leslie, P. H. (1945). The use of matrices in certain population mathematics. Biometrica, 33, 183. Ravenholt, R. T. and John Chao (1974). World fertility trends. Population Report, Family Planning Programs, Series J, No. 2. Song, Jian, and Jingyuan Y u (1981). On stability theory of population systems and critical fertility rates of women. Acta Automatica Sinica, 7, 1 (in Chinese). Wang. Huanchen (1980). Population state and population dynamics, Journal o f Xi'an Jiaotong University, 14, I (in Chinese). Zhu, Yue-xin (1981). Research on the parameter identification of a population model of China. Master's degree thesis, Xi'an Jiaotong University, China (in Chinese).
n-I
Hi,,(k)w(k + n - i - 1) = H~(k}w,(k).
(55)
i=O
Substituting (54) and (55) into (53), we then obtain (12a}. 2. Noting that yo(k + j ) = b'(k + j)xlk + j )
= b'(k +j){Hj(k)x(k) + H~(k)yoj(k) + H~(k)w~(k)}, j = 1, 2. . . . . n and b'(k +j)H~(k)yoi(k) = O, j = 1, 2. . . . . n < 0 a n d from (15) and (22a, b), we have H~(k)wl(k ) = H~a(k)wx(k),
(56)
j = 1, 2 . . . . . n.
Hence yo(k + j ) = b'(k + j)Hj(k)x(k) + b'(k + j ) H ~ (k)w,(k), j=l,
(57)
2,...,n.
Substituting it into (14), we have y6,(k) = [b'(k
I
I
I
I
'
1)H,_ t(k)x(k)l...~ ' b'tk)Ho(k)x(k)lOl.., i0]
+ n -
I
I
+ [h'lk + n - I)/-P'~ (k)w,(k)[...l
t
b'lk + 1)H~t(k)w,,(k) ~)0i...i 0]' I
I
I
j
•
I
+n-l),...i
+
w',(k) [ / ~
t'(k)b(k
+
n
-
Io]
1 )),' '" II
H~i(k)~(k + 1)i01. ..lO] ' A P P E N D I X A. A DERIVATION OF (12) A N D (13)
x'(k)B',,(k) + w',(k)B~'(k)
I. Letting u(k) = qo,,(k)yo(k) + w(k)
(51)
the transpose of which is just (13a)' 3. From (18), we have
from (9) we have
2,
x(k + 1) = H(k)x(k) + u(k)
sin(k)
x(k + 2) = H(k + l)H(k)x(k) + H(k + l)n(k) + n(k + 1)
=
,,,,k+2,.
__ ,,,,,,k.
L;(k+ n)
Li'x(k + n)J Li'H.(k)J
l'Ht(k)yo,(k)
tl-I
x(k + h) f Hh(k)x(k) + ~ Ho,(k)u(k + h - i -
x(k)+
=
Fl'H~(k)w,(k)]
1),
i=O
h = 0. I ..... n
(52)
where Hh(k), Hi.h(k) are given in (20) and
i'H~(k)yo,,(k).] Li'H•(k)w.(k)J Now, an arbitrary component of the second vector on the righthand side of (58) can be written as
-1
y [.] = o.
l'H~(k)yoh(k ) = q_~(k + h - l)yo(k + h - 1) + ...
i:0
+ q-I(k)qo(k + 1)...qh-2(k + h - l)yo(k)
Letting h = n and combining (51) and (52), gives
= l'H~l(k)yo,(k),
x(k + n) = H.(k)x(k) n-I
where H~l(k) is shown in (27b). Putting (59) and (56) into the second and third vector on the right-hand side of (58), respectively, we then obtain
+ ~ Hi.,,(k)qoa(k + n -- i -- 1)yo(k + n - i - 1) i=O ¢1-1
+ ~ H~.,(k)w(k+n-i-
1)'
(53)
i=O
s,(k) = F~(k)x(k) + G,(k)yo,(k) + L~(k)w,(k)
H~.,(k)qo,,(k + n - i - 1)yo(k + n - i - 1) =
]
Fd(k)
1
[ 0.p I I /I i=~'tr//(k + n - i +j)yo(k -~ n - i - 1)1~0 i . . . 1,~7 0 /j, _ i
i + t ..... m
Ld(k+n-l)J
and therefore + H~,.(k)qo.(k + n - i - l)yo(k + n - i - I ) = H~,(k)yo.(k).
F
(54)
L/~'(k+n-llx(k+n-1
{I -,lo,,(k)]yo(k) {! - r/o,(k + 1)}yo(k + 1)
LIt -
/=o
Noting (15)and (22a), we have
(60)
where Fs(k), G,(k) and L,(k) are (26a, b, c), respectively. Similarly, we have
From (2), (20), (3), and then (14) and (21), we have
0 ..... i - - 1
(59)
h = 1, 2 . . . . . n
l
,J (6l)
rlootk ÷ n - 1)}yo(k + n - 1)
where the vector~(k + j ) , j = 0, 1..... n - 1, is in (27a). The first vector on the right-hand side of (6l) represents the total death
The p a r a m e t e r identification of a population model of China population in state-transition process; and the second one the death number of babies born in the very year. Putting (12a), (27c) and (14) into (61), we have
d.(k) = / { (k + L'~'(k + n - 1)H._dk
+
because, at ~ ffi 0 u(t - ~)~.o = u(t) = p(O,t).
Hence the function
)]
x(~,t) =
po(~,t)d~,
~¢[-1,m-
x(a, k) = xo(k),
(67)
From (63)-(67), the function defined as ~¢[-1, m-1]
tp(~,t)=x(~+l,t+l)/x(~,t)
Substituting (591 and (56) into the above, we then obtain d.lk) = F~(k)x(k) + Galk)yo,,(k) + L,~(k)w,,(k)
(62)
where F.(k), Gd(k) and L,~(k) are in (26b), (26d) and (26f), respectively. Substituting (12a), (60) and (62) into (17), and noting (25), we have
/ = I r'(k) I x(k) ÷ l~'(k)lY°"(#)
Ld.(k) J L~',(k)_l
L~Ak)J
l
w.(k)l LL,(k) J |L,(k)
ffi Co(k)x(k) + D.(k)Yo,,(k) + Ro(k)w.(k)
which is just (12b). Equation (13h) is obvious. APPENDIX B. A PROOF OF THE PROPOSITION We assume it to be reasonable that continuous mathematical models are used to describe population process~. This is true, i n particular, when the total number of population individuals is large e n o u g h - - that is the way of Chinese population. In continuous population models (Focrster, 1959; Song and Yu, 1981), the state variable function is a continuous distribution function, with respect to
p(~,t)>o, ~ C0, m] which represents the population number per unit age at age ~ and at moment t; and the input variable function is a continuous distribution function, with respect to t
u(t)>O,
t~(-~,
+oo),
which represents the birth-population number per unit time at moment t, and, when ~ ffi 0
u(t) = p(O, t).
~
is just the age-specific survival function ~(~) at moment t we are looking for. This is because, except for x(~,t) = 0, the function tp(~, t) is obviously differentiable with respect to ~; and, for V~ ¢ [ - 1, m - 1 ], the function ,p(~, t), in the sense ofthe discrete model is just the survival at age ~ and at moment k; and, in particular, ~p(~,t)k=,=q.(k),
a = - l , O, .... m - 1 .
tlk
APPENDIX C. ON DETERMINABILITY OF PARAMETER ESTIMATES By means of the least-square method, the identification of both the age-process and the birth-process is transformed into a minimization problem. Obviously, especially under the condition of those very limited observed data, it is essentially important that whether the solution of the minimization problem, i.e. the parameters under estimation, be unique or not. If not, the process will not be identifiable anymore by those data. Conversely, if the solution is unique, it will be the most probable solution of the parameters determined by the observed data. Although some other arbitrariness still exists, we call this process determinable under the observed data, and regard this solution as the final result. 1. Definition. For certain observed data available, a process is said to be determinable under the observed data if its model parameters can be uniquely determined by those data. 2. Determinability o / t h e age-process. Determinability of the age-process under the data structure here depends on whether the optimal solution of (31) with respect to ~ is unique or not. However, we have Theorem 1. The optimal solution of the nonlinear optimization problem (31) with respect to ~ is unique. Proof. It is enough to prove that, without the statistical data i, a and ~,(k), this solution of (31) would still be unique. In fact, letting ~J.(~) 8q
we can then analYtically find the unique optimal solution qo as a-I
7° = ~,+1(k + n)/~o(k + n - a - 2 ) / r ]
which is just yo(k) in (3b) when above equation, we have
:o
qJ,
j=-I
t+ 1
yo(t) ffi
(68)
~o ffi :~o(k + n ) / ~ o ( k + n + 1)
Under this condition, the aggregate number of the birthpopulation born in [t,t + 1) can be written as
yo(t) =
(66)
a = O, 1. . . . . m
x ( - I, k) = yo(k).
°
+
(65)
and, from (63)
~ ' ( k + , - 1)n:-'Ik)yo.,- ~(k)J
z,.(k + n) |s.(k)
1]
is differentiable with respect to ~ and is the aggregate population number in [~, ~ + 1) at t. When ~ ffi a and t ffi k, then
~.'(k + 1)Hl(k)yodk)
I
427
a--0,1 ..... n-2
u(~)d~ a-1
t ffi k. L e t t i n g
~
--- t - ~ in the
~°ffi~.+1(k4n)l~.-.+,(k)/
H j-e--n+
~, 1
a=n-l,n-2~
u(t - Od~.
(63)
-1
We now define an auxiliary function
~
~E I-l,0)
(64)
which is obviously continuous with respect to ~ ¢ [ - 1 , m ]
.... m - l .
The detail of the proof is omitted here. 3. Determinability o f the birth-process. If the parameter vector h is taken as a direct object of the parameter estimation, the birth-process described in (33) will obviously be not determinable; but we have Theorem 2. The birth-process described in (34) is determinable if and only if rank X ffi ~ - 0 + I.
428
YUE-XIN Z r t c a n d B A b W U W A N
Proof. If
birth-process is rank X ~ ( - 0 +
1
rankX<~-0+
1
rank X ~> r.
there must be
Proof. Discretizing
rankX=~-0+l there must be either
v) b(O +
I, v), .... b((, v)].
Since b(~, v) is continuous and differentiable in expressed as
v-space, it can
be
(69)
b(v)=b(vo)+{~}'Av+O(Av)
(71)
I
rank [X II~'] = rank X or
b(~, v), we have
h'(v) = [b(0,
and (34) must have infinitely many solutions of b and is not determinable. Conversely, if
(70)
I
wherev : Vo + Av; 0(Av) is higher order infinitesimal vector than Av. Neglecting 0(Av) in (71), and substituting it into (34), we have
AAv
rank [X ll~ ] > rank X In the former case, (34) is consistent and has a unique solution of b. In the latter ease, (34) is contradictory and its regression solution b = (X'X)-
= g
(72)
where • f~b'(vo))
A= A ~ ;
'
(73)
tX'i,
exists uniquely because of (69). /..emma 1. The birth-process is not determinable if I < ( - 0 +1. /,emma 2. The birth-process of a stationary population under this data structure is not determinable. These two lernmas both are obvious because
g = ~ -
Xb(vo)
with dimension I × r and I × 1, respectively. For the vector Vo of arbitrary choice, the matrix A and the vector g can both be regarded as known, and the linear regression problem (34) with respect to b has been transformed into (72), with respect to Av. From Theorem 2, the neoessary and sufficient condition of the determinability of (72) is rank X = r. From the theorem of rank of a product matrix, this theorem is then proved. Remarks. 1. From the above proved, it follows that
( r a n k X)max = l < ~ - 0 + l
rank
and rankX = I respectively. 4. If b(~, v) is the relation among components of the timeinvariant vector b where v is the independent parameter vector with dimension r; and b(~, v) is continuous and differentiable in vspace; then we have:
Theorem 3.
A necessary condition of determinabifity of the
= r
(74)
is another necessary condition of determinability of the birthprocess. Because the condition (74) is often satisfied easily here, the condition (70) also has some suff~ientness, or rather, the determinability of the birth-process is mainly determined by (70). 2. When the relation among components of the time-varying parameter vector b(k) is expressed as (46), the conclusion about the determinabifity of the birth-process is the same as (70), and its proof is not presented here.
Automatica. Vol. 20, NO. 4. pp, 415-428. 1984 Printed in Great Brilain.
The Parameter Identification of a Population Model of China* YUE-XIN Z H U t and BAI-WU WAN++
A mathematical form for the survival function drastically reduces the degrees of freedom in the existing population model and makes the solution of this parameter identification problem, under very limited statistical data, practicable. Key Words--ldentification; least-squares estimation; modelling; optimization; parameter estimation; population model.
death population in several successive years to estimate the age-specific mortalities directly (Boyarsky and Shusherin, 1955; Keyfitz, 1968). But these data do not exist in the ordinary census in China. That is why population model parameters for China have not yet been obtained. Hence, the study of estimating the population model parameters on the basis of the statistical data currently available is important. Our goal is, on the basis of the statistical data available, which include the data from censuses in the past 25 years, to study the parameter estimation problem for a population model of China by system identification, and to try to acciuire some useful model parameters. By comparison with common physical systems, one of the features of population model parameter estimation is that past statistical data, no matter how limited, are the only window from which the behaviour of the population process can be observed. The statistical data from the Chinese population system which are available for parameter estimation are very limited. Another feature is that the data-sampling time-periods of different output data of the population system are not identical, so the model structure assumes a very complicated form. Thirdly, the population model is a bilinear control system model with a rather high dimensionality and quite a lot of parameters to be estimated. Finally, the irreversibility of the population process means that the population forecast and, therefore, the prediction of model parameters themselves, have great importance. All these features give rise to the difficulties in the parameter estimation of a population model of China. In this paper five steps are given: (1) we choose a discrete mathematical model to describe the population process; (2) outline the available statistical data and their preprocessing; (3) transform the common model structure into a form
A l ~ ' a e t m A mathematical model for predicting and controlling the future trends of the population of China is given. The model used is a discrete bilinear system with 101 dimensions. The model is decomposed, transformed and its parameters are estimated by nonlinear optimization and curve-fitting techniques. Results for both time-invariant and time-varying parameter estimation are obtiined. Finally, a successful validation of the population prediction over the past 25 years is given, in which the predicting accuracy ofthe total population for each year is within about 1 ~. This is despite the fact that the available statistical data for these years, which include two censuses, are limited, incomplete and contain numerous contradictions. INTRODUCTION
CONTROL theory has been finding important applications in the study of the population process. With the deepening of these quantitative studies, one problem in particular has arisen--the parameter estimation problem of a population model, because it directly affects the credibility of the process forecast, the reliability of its synthesis and, therefore, population policy decisions. However, the population model parameters, i.e. the age-specific survivals and the age-specific fertility rates, have all traditionally been evaluated by means of demographic methods which require much more detailed statistical data than ordinary census data. For example, those demographic methods need the statistical data for the age-specific distribution of
* Received 30 March 1982; revised 2 December 1982; revised 13 October 1983. The original version of this paper was presented at the 6th IFAC Symposium on Identification and System Parameter Estimation which was held in Washington, D.C. during June 1982. The published proceedings of this IFAC meeting may be ordered from Pergamon Press Ltd, Hcadington Hill Hall, Oxford OX3 0BW, U.K. This paper was recommended for publication in revised form by associate editor J. Mendel under the direction of editor H. Kwakernaak. t Department of Electric Engineering, Xi'an Institute of Optics and Precision Mechanics of Academia Sinica, Xi'an, Shaanxi, People's Republic of China. :~Institute of Systems Engineering, Xi'an Jiaotong University, Xi'an, Shaanxi, People's Republic of China. 415
416
YUE-×IN ZHU and BAI-WU WAN
corresponding to the real data structure; (4) decompose it further into two subsystems--the age-process and t he birt h-process, and identify them separately to simplify the problem; and (5) formulate the parameter estimation problem and convert it into an optimization problem with dimension 101. In spite of this, it is still too difficult to estimate the model parameters of either the age-process or the birth-process directly because the statistical data are too limited. So we find some relations among the components of the parameter vectors under estimation by a curve-fitting technique, and thus simplify the nonlinear optimization problem by reducing its dimensionality. We finally show that the optimal time-varying parameter estimates, as well as the time-invariant ones, can both be found, and, give a satisfactory validation. MODEL STRUCTURE AND ITS TRANSFORMATION
Choice o f the mathematical model
Several kinds of mathematical model have been developed to describe the population process. They can be briefly summed up as follows: first, continuous mathematical models with continuoustime, continuous-age and continuously-distributed state variable, which are expressed as partial differential equations (Foerster, 1959; Song and Yu, 198 l). Second, discrete models with discrete-time and discrete-age variables and lumped-state variables, which are mainly expressed as a Leslie matrix or one of its variants (Leslie, 1945; Kwakernaak, 1977; Wang, 1980). Third, a model characterized by a finite-state continuous-time semi-Markov process (AIj and Haurie, 1980). Given a large scale of China's population system and its statistical data structure, we prefer a discrete model to describe the population process. The general structure of this model is (Kwakernaak. 1977; Wang, 1980) x(k + 1) = H ( k ) x ( k ) + qo~(k)B(k)x(k) + w(k) (1) where k is a discrete-time variable, in year, with origin at 0AD (the moment k, refers to the appointed data-statistical time-moment in the chosen year kAD)
a discrete-age variable. The oldest age (considered in the model) is m = 100
H(k) =
.
i]
~7~(k) ....
(2)
r/,._ l(k)
is a Leslie matrix with dimension (m + l )2, in which all the elements are equal to zero except t/a(k), a = 0, 1 ..... m - 1, on the first lower-left subdiagonal of it. In (2) G ( k ) = xa ÷ i (k + 1 )/x,(k),
a = O, 1. . . . . m - 1
(3a) ~loa(k) = xo(k + l)/'yo(k) = r/_l(k)
(3b)
are called the survival at age a and the baby survival in the very birth-year (Wang, 1980), respectively, at the moment k, in which yo(k) is the aggregate number of babies born in [k, k + 1). For simplicity, we define t/'(k) = [t/oa(k) r/o(k)..... r/,._l(k)]
(4)
as the age-specific survival vector with dimension m + 1. And B(k) = [b(k) i O(m+l)×m]'
is the birth-matrix where b'(k) = [bo(k)
b l ( k ) ..... b,.(k)]
(5)
is the age-specific fertility-rate vector and ~0, O > a o r a > L bo(k) = [yo.Ak~/x~(k), 0 < a < ~
(6)
where yo.~(k) represents the number of babies born by x~(k) in [k,k + I); and 0 and ~ are the youngest and the oldest fertile-ages, respectively. In this paper, 0 = 17, ~ = 50. The vectors q(k) and b(k) are the model parameters to be estimated. The term b(k) is also usually written as b(k) = k:(k)fl(k);.(k)
(7)
where k : ( k ), fi{k) and ).(k) are the ratio-coefficient of women in fertile-population, the average fertilityrates of women and the fertility-pattern vector, respectively"; and
x'(k) = [xo(k)xt(k) . . . . . .'¢~(k)] (8a)
fl(k) = V h,,(k),k/.(k}
and
a =0
2'(k) = [2o(k)2~(k) ..... 2~(k)l
w'(k) = [wolk)wt(k) . . . . . w.,(k)] are, respectively, the population state vector and the population state disturbance vector at the moment k, the prime indicates a vector transposition; x~lk) and wAk), a - 0, l ..... m, are the aggregate number of people and the number of disturbance people, respectively, in age-region ,~a.a + I ). The a denotes
2,,(k) = b,Ak)/{kr(k)fl(k)~.,
(8b)
a = O. 1. . . . . m
(8c) and m
Y :.,,(k)= 1.
The parameter identification of a population model of China Obviously, the birth-process can be more clearly controlled and described, and therefore identified in terms of it(k) and fl(k) than in terms of b(k).
417
z',r(k) = [yo(k)yo(k + 1). . . . . yo(k + n ,
1) 3 (16)
and [
z~o(k + n) = [x'(k+ n) , s;(k) ',,d'.(k)] Transformation of the model structure Equation (1) is not suitable for the parameter estimation because it is difficult to identify this bilinear system directly; its discrete-time period is different from the real data-sampling time-periods; and there is no output equation corresponding to real statistical data. Hence, the model structure (1) must be transformed. We first divide the model (1) into two subsystems: the age-process subsystem x(k + 1) = H(k)x(k) + r/oQ(k)yo(k) + w(k) (9) and the birth-process subsystem yo(k) = B(k)x(k)
(10)
y~(k) = [yo(k) 0,..., 0 ]
(11)
where
is the birth-vector with dimension 1 x (m + 1). These two subsystems are common linear systems, directly coupled with each other, and both are concerned with only one of the two parameter vectors ~/(k)and b(k), respectively. So we can identify them both separately to avoid identifying the bilinear system directly (Eykhoff, 1974; Isermann, 1979). Next, assuming the data-sampling time-period about the population states to be n years, we derive the model structure corresponding to the data structure, including the output equations, as follows:
are the output vectors of the birth-process and the age-process, with dimensions 1 x n and 1 x (m + 2n + 1), respectively, s',(k) - [s(k + 1)s(k + 2). . . . . s(k + n)]
[d(k)d(k + l) ..... d(k + n - 1)] (19)
d,,(k)
are the total population vector and the total death population vector in [k, k + n), respectively, where s(k + j) and d(k + j - 1),j = 1, 2 . . . . . n, are the total population at k + j and the total death population in [k + j - 1, k + j), respectively. The parameter matrices in (12) and (13), can be briefly written as h
Hha,(k) = l-I H(k + p - j ) ,
(20a) in particular
Hp,r(k ) = Hp(k),
H, ( k ) = H(k),
where I,"+, is the identity matrix with dimension (m + 1)2; and
Hi(k) = dia
--ri
L.i=
th(k + h + i)
- 1
0 0
17[ th(k + h + i - 1)..... i= -1
1
(12b)
for the age-process subsystem, and yo,(k) = B,(k)x(k) + B~,(k)w,,(k)
(13a)
z~r(k) = Eyo,(k)
(13b)
h_2 I-I rh(k + i +
1
1) 0 , . . . , 0 , h = I, 2 . . . . ,n
i=-1
h-
1
h,...,m
for the birth-process subsystem where
with dimension (m + 1)2; and
y~,(k) = [yo(k + n - 1)yo(k + n - 2).....
H~(k) = [Ho.h(k) ,[ H, .h(k)', .....
(21)
' [Hha.h(k)], h = 1, 2 . . . . . n
(14)
(22a)
and
and w'.(k) = [w'(k + n - 1) ,l w ' ( k + n - 2 ) ]
(20b)
Ho.p(k) = Ho(k)= 1,,+1
z,,,,(k + n) = Co(k)x(k) + Do(k)yo,,(k)
yo(k)O ..... 0]
h = 1,2,...,p ~< n
j=l
(12a)
+ Ro(k)w,,(k)
(18)
and
x(k + n) = H,,(k)x(k) + H~(k)yo.(k)
+ H~,(k)w,,(k)
(17)
! h H~l(k) = [0(,"÷1)x.,-h),Hw(k)]. h = 1,2 . . . . . n
, ..... w'(k)]
(15)
are the total birth-vector and the total disturbance vector in [k, k + n), with dimensions 1 x (m + l) and 1 x n(m + 1). respectively
(22b) with dimensions (m + l) x h(m + l) and (m + 1) x n(m + 1). respectively, and, when h = n
H~, (k) = H~(k)
(22c)
418
YUE-X~N ZHU and BAz-WU WAN I
B'.,(k) = [H~,_t(k)b(k + n - l) I
/{oo(k) = diag[l - qoo(k + n - 1) t
t
1 - qoo(k + n - 2) .....
H/,_2(k)b(k + n - 2)~r.... i I !
I ))"
H;(k)b(k))O~
..,i 0]
1 - qoo(k) 0 . . . . . 0]
(23)
and
(27c)
both with dimension (m + I)2; and
B~,(k) = [ H ~ ( ~ ( k ) b ( k + n - 1) . . . . . i
H~lt(k)b(k + 1)
0
r=
i . . . . . 'n, O ]
(24)
with dimensions (m + 1)2 and n(m + 1) x (m + 1), respectively. The observing matrices are (25a)
C'.(k) = [H,,(k) l f .(k) I FS(k)
[H,(k)~G~(k)~Gi(k)]
D'o(k) =
(25b)
both with dimension (m + 1) x (m + 2n + 1), and R'o(k) = [I-t"~ (k) l L,(k) l L'd(k) ]
(25c)
and
L1' ,1',, [o 0,0,×,,,_.,-1) 1
E =
(25d)
with dimensions n(m + 1) x (m + 2n + 1) and n x (m + 1), respectively. In the above, iI , l F'~(k) = [ n ' , ( k ) l Hz(k)l ~). . . . ..' ttH;(k)l ] d
(26a) i
F'n(k) = [H'o(k)l~(k) H~(k)l~(k + l)! ..... i H~_l(k)l~(k + n - 1)]
(26b)
G'~(k) = [H,~l(k)l ,')H2t(k)l iI , .... iIH,. l(k)l ]
(26c)
G'.,(k) = [0 H ~ ( k ) l ~ ( k + 1), .... , H~L-l(k)l~(k + n - 1)] + I~o~(k)E'
(26d)
all with dimension (m + 1) x n, and I
,2
I
I
t'~(k) = [H~ll(k)l ~H~t(k)l I ..... ~H~,(k)l]
(26e)
..,i i
I ~L L'~(k) = [O[Hwl(k)~(k + 1 ) II, .
,
' " - ' (k)l~(k - + n Hw~
1)]
(26f)
both with dimension n(m + 1) x n. Also in the above #'(k + j ) = [1 - qo(k + j ) 1 - r/t(k + j ) .....
j = 0 , 1 ..... n - 1 (27a)
l-qm_~(k+j)l],
with dimension 1 x (m + 1); and
[1
1,...,1]
the dimension of which corresponds to that of a matrix multiplied by the vector, 1, or that of another vector multiplying the vector 1 in equations. Equations (12) and (13) are typical and canonical linear control systems but their internal structure is rather complicated. Their derivation is given briefly in Appendix A. Matrices H,(k), H~(k), H~(k), Co(k), Do(k) and Ra(k) in (12) contain only t/, while matrices B,(k) and B~,(k) in (13) contain both q and b. So the age-process (12) can be identified solely by the parameter vector q alone. Hence, we can identify the birth-process after having identified the ageprocess. TRANSFORMING PARAMETER ESTIMATION PROBLEMS INTO OPTIMIZATION PROBLEMS Available statistical data and their preprocessing The original statistical data available for identification are not only limited and incomplete but also contain numerous contradictions. Therefore, we had to preprocess them by: standardizing, their data-collecting time-moment; converting the regional data into the national ones; converting and fitting the baby mortality data; investigating the precision of all data used; and further processing the parts of those data containing obvious contradictions and serious errors, and so on. Then we obtained all the preprocessed data, and their precision distribution, which consist of: (a) ilk), k = 1953, 1964, 1975, where the symbol " denotes statistics of the vector under it, of which i(1953) and i(1964) come from the 1953 and 1964 national census data, respectively; while i(1975) comes from the regional census data of 1975 and have been converted into national ones. Thus, there are only three sets of sampling data of the population states at three different moments in two data-sampling time-periods, both with n = 11 years. (b) g = [stk)s(k + 1)..... s(k + 25)] a ' = f l l k ) d ( k + 1)..... d(k + 24)] ~' = [Yolk)Yolk + 1). . . . . yolk + 24)]
(28)
-1
~~0, 1-I r/ilk + h + i) l ..... i
H i t ( k ) = diag E0
t=--
0..... n-h-1 n-2
i
n -
t
y~ ~,(k + i + 1)~0~ . . . . . " ',"
i=-1 n -
and
I
n
1
,...,
I
m
h
...
where k = 1953. These all come from routine annual statistical data, and their data-sampling time-period is one year. The formulation o f the parameter estimation problem
h = l ..... n (27b)
From (17), we see that the state vector x(k + n) is contained within the output vector z.o(k + n) of the age-process (12), which has been arranged on
The parameter identification of a population model of China purpose, and therefore, the state equation (12a) is contained within its output equation (12b), too. Hence it ~s enough for the age-process (12) to be studied through only its output equation (12b). Considering the recurrence properties of discrete models, we can represent the age-process as the model with input vectors x(k), yo.(k) and w.(k), and the output vector z..(k + n), shown in the dotdashed block in Fig. l(a), if there were no symbols ^ and - above the vectors. As for the birth-process (13), it is really a set of algebraic equations. Because of the special structure of matrix E, there is essentially no difference between the two vector, yo.(k) and z.v(k ). Hence, the birth-process can also be represented as a model with input vectors, x(k) and w.(k), and output vector yo.(k), which is shown in the dot-dashed block in Fig. l(b). Now we compare the mathematical model with its own real process for the age-process and the birth-process, respectively. As shown in Fig. l(a) and (b), respectively, we take the input vectors of a real process corrupted by measurement noises (i.e. statistical errors) as the input vectors of its model, and then use the output response of the model and the output vector of the real process corrupted by measurement noises to form the error vector as
er(k) = ~o,,(k) - ~o,(k) = B,(k)~(k) + B~,(k)f%(k) - ~o,(k)
Transforming the parameter estimation problems into optimization problems In the case of the age-process, considering the age-specific survival vector # as time-invariant during every time-period, and in the sense of weighted least squares, we transform the parameter estimation problem of the age-process into a nonlinear optimization problem, for every timeperiod respectively, as follows:
Co(k)f~(k) + D,,(k)~o,,(k)
(29)
+ Ro(k)~k,,(k) - ~,,°(k + n)
with dimension ( m + 2 n + l ) x 1 for the ageprocess, as shown in Fig. 1(a), or as
vz(k+n) w n(k) x(k)
[
Yon(k)
/
,.(,)
I
Zno( k +n J
,+++Z~lk+n)
't '1"" Ji;rm(k + n)
-leo(k+n)
ReoL age -process
,
I
Colk)
Do(k) Math. model of the age prooess
wn(k)
X(~)
,
'* '1 Ro( , •
]
[
ReaL
(a') _J
~.(k}
,.,,,
II I
L
i
vy(k)
~(~+ ;o. (k)
.
birth - process
;r =[ .:(k, I+.°o,+ v,lk),~÷ x(k)
Yon(k)
~l-*y (k)
(b) Math. rn~et of birth-process ----.-J
FIG. 1. The formulation of identification problems: (a) for the agc-proo:ss, (b)for the birth-process.
AUTO20:6-'C
(30)
with dimention (m + l) x 1 for the birth-process, as shown in Fig. l(b). In the above two equations and in Fig. 1, the overbar denotes a response of a model to its inputs. The vectors v=(k), vy(k), vw(k) and v.(k + n) are the corresponding measurement noises. Note that the error vectors, (29) and (30), are both expressed as functions of the parameter vectors r/ and b under estimation, respectively. What either of the two parameter vectors is taken as determines the corresponding error vector and, so determines the closeness between the model and its real process and, therefore, the effectiveness of the model, and vice versa. Hence, the specific formulation of this parameter estimation problem here ought to be: choose the model parameter vector so that the corresponding error vector, in a sense, is a minimum. It is worth noticing that these two parameter estimation problems are both problems in which the input a.nd output-signals are corrupted by measurement noises.
eo(k + n) = ~,(k + n) - $,,(k + n) =
419
420
YUE-XIN ZHu and BAI-WU WAN b'(k + j ) = [bo(k +j)bo+,(k +j) ..... bdk + j ) ]
,.a~-,
+ 2n + 1 e.(k +
(31)
where R" + 1 is a real vector space with dimension m + 1; J~(q) the objective function with the same dimension as the mean square error; e,(k + n, q)just the error vector e.(k + n) which is a function of q; Q the weighting diagonal matrix determined by the precision distribution of corresponding data. Thus the parameter estimation problem of the age-process has been transformed into a feasible but rather complicated and intrinsically nonlinear optimization problem with dimension 101. As for the birth-process, since b(k) appears linearly in B,,(k) and B~,(k), and also in the error vector ey(k), as seen from (23), (24) and (30), where q(k) has been considered as known in identifying the birth-process as mentioned above. So the parameter estimation problem of the birth-process is really the unique regression-solution problem of the following linear matrix equation (Eykhoff, 1974)
$o.(k) = B.(k, b),t(k) + B~(k, b)#.(k)
(32)
with respect to b, where B,.(k, b) and B~(k, b) are just B,,(k) and B~(k), respectively, and are linear functions of b(k). Since q is known, (32) can be written as
~o,(k) = [x'(k + n - l)b(k + n - 1},...I
01 where x(k + j ) , j = 1,2 ..... n - 1, can be derived from the real statistical vectors :~(k), w.(k) and the known vector q according to (12a). Considering the limitation of the birth-statistical data, we identify the birth-process not only within a s t a t e - - d a t a sampling time-period but all over the 25 years for which we have birth-statistical data, as shown in (28). And, for simplicity, when treating the birthprocess, we cut out the redundant zero elements in vectors yo.(k), b(k) and A(k), and represent them with the same symbols as above. Thus, the parameter estimation problem of the birth-process can be finally transformed into the unique regression-solution problem of the linear matrix equation, either ~ ' = [x;(k + l -
l)b(k + l -
and x;,(k + j ) = [xo(k +j) Xo+,(k +j) .....
x;(k+j)],
j = 0 , 1 ..... l -
both with dimension 1 x ( ( - 0 matrix X ' = [xb(k + l -
1
(35)
+ 1); and the
l)!xb(k + l - 2) i... I Xb(k)]) I (36)
with dimension (~ - 0 + 1) x I. In the same sense as the above, xb (k + j) and therefore X are considered as known here. IDENTIFICATION OF THE AGE-PROCESS
The solution of the optimization problem (31) and difficulty with the parameter estimation In order to solve the optimization problem (31), we used a TQ-6 computer (a large computer made in China); and used the Davidon-Fletcher-Powell optimizing algorithm which requires an analytical expression for the gradient of the objective function. We solved the optimization problem for the two time-periods, i.e. 1953-1964 and 1964-1975, respectively, and then obtained two optimal estimates of the age-specific survival vector. A typical example of these two results, for 1953-1964, is shown by the discrete-point set in Fig. 2. The other results are similar. Obviously, such a result is impractical. In other words, these optimal parameter estimates, from a direct solution of (31), are too dispersive either to be taken as the final result of the parameter estimation, or to be further processed by means of the usual techniques of graduation of data. This is just the difficulty with this parameter estimation by a direct solution of optimization problem (31). Obviously, in view of system identification, that is because of too limited statistical data.
Analytic relation among the age.specific survivals Concept. In the above as well as usual discrete population models, all the components of the vector ~/have been regarded as independent. However, they are interrelated through the human-fife-process in a
l)i.., i i¢;(k)b(k)] i
°
oI
(33)
:
with respect to b(k + j),j = 0, l . . . . . l - 1, which are assumed to be time-varying; or 9
= Xb
(34)
"~. oaok I
with respect to b which is assumed to be timeinvariant. In (33) and (34), k = 1953; l = 25 denotes the number of the birth-statistical data used ; ~ is just the birth-statistical data vector (28), with dimension lxl
" ""
070~-
.
• \
i
~°6° I-F-
:~_ 050
0
.,~
~
r
I
I
I
1
I
I
t
tO
20
30
40
50
60
?0
BO
90
Age
(yeors)
FIG. 2. A direct solution of (31).
,
The parameter identification of a population model of China society so that among them there would exist a relation which is the very thing we are interested in. On the other hand, strictly speaking, the survival in this paper as well as in demography is, in nature, a probability of survival. However, the surxlval defined in (3), i.e. the survival quotient for individuals, is only a frequency of occurrence of survival under the condition of finite individuals. Obviously, the probability of survival is determined by social factors, medical conditions and so on, while the closeness of the frequency to its probability is determined by the number of individuals under discussion. The relation among the age-specific survivals discussed here is just the dependent relation among the individual probabilities of survival for the different age-set population. Hence, this relation is deterministic. For simplicity, we will still use (3) as the definition of the age-specific survivals, but understand it in the same way as the above.
The survival function. On the basis of continuous population models existing, we now have a proposition as follows. Proposition. For a set of the age-specific survivals in a discrete model at a moment k, r/o(k), a = - 1, 0 ..... m - 1 , there must exist a continuous and differentiable function r/(~), with respect to the continuous age-variable ~, which, in the sense ofthis discrete model, is just equal to the age-specific survival at the age ~ e [ - 1, m - 1 ] and at the same moment k. In particular, when ¢ = a, the function r/(¢) satisfies r/(a) = ~/,(k), a = - 1 , 0 . . . . . m - 1.
421
and invasion of diseases, reach some equilibrium, and the survival reaches a maximum. After that, the increase in the factors resulting in death exceeds the influence of the increase in adaptability, and the survival decreases slowly over a rather wide agerange. After about 50 years of age, for example, senescence begins to increase, and after some age, for example 60, the survival starts decreasing quickly. The greater the age, the faster the decrease of the survival until the survival vanishes. That is the rough shape of the survival function, and therefore, its degree of freedom has been restrained.
An analytic relation. From the above, the teenager section of the survival function can be approximated by an exponential function; the middle-age section by a linear function; and the old-age section by a composite power and exponential function. Summing these, we then obtained an approximate expression for the survival function with seven independent parameters
-
-
ctaexp(--at2~ ) +
expt - 0,• I ~e [-1, m)
) (38)
where
is the independent parameter vector with dimension 7; the factor in the linear term
(37)
The proof of the proposition is given in Appendix B. Obviously, the function ~/(¢) is just the relation we are looking for, and we call it the age-specific survival function, or simply, the survival function. In order to find the specific expression of this function, let us look at the real statistical law of the age-specific survivals.
Statistical laws. Many real statistical results show that the age-distribution of the survival has the following laws: The survival of a group of new-born babies is generally lower because their existing state changes rapidly and they are poorly adapted to their new surroundings. As their age increases, their adaptability increases quickly and so does their survival. However, the older they are, the lower the rate of increase of their survival becomes. On the other hand, with the increasing age, the probability that diseases and other external harmful factors will affect people rises slowly but surely, and almost linearly. At some age in the teens, the two factors, the increase of adaptability to external surroundings
J is chosen so that it does not greatly influence the value of ~(~, e) when ~ < 70, but it gives ~(~, ~) the following analytic features
> 0,
[ - 1 , oo)
and as ~ --* oo, r/(~, at) --. 0 monotonically.
Effect. Is it effective that the survival function (38) is used to express the age-distribution of the survival? Obviously, this can not be proved theoretically and can only be verified by practical survival data. An effective survival function should be able to approximate various kinds of age-specific survivals closely enough, by means of changing the independent parameters in it. In order to verify the function (38), we used it to fit several kinds of agespecific survival data which came from real, published statistical results and had been graduated. This is a curve-fitting problem, i.e. the solution of the nonlinear optimization problem
422
YUE-XIN ZHU and BAI-WUWAN min J,(,,) =
where R 7
/ [ m - - ~ ],q(at) - @H~1
min Jo{q(=))
(39)
~t~R 7
,
is real vector space with dimension 7
~'(,')= [ q ( - l , ~ )
q(O,~)...q(m- 1,=)]
(40)
is the fitting function vector; and
#'= [qo,, @o.-.0,.-i]
(41)
is the age-specific survival data to be fitted as a vector. This was done with the DFP-algorithm on a DJS-130 computer (a small computer made in China). The results of fitting this function to many sets of real statistical data show that the survival function (38) is satisfactory. As an example, Fig. 3 shows a result of fitting (38) to a kind of the age-specific survival data, which are published in (Song and Yu, 1981). The largest fitting errors are within I%0, when age <57; and within 1 ~ , when age <70. And the fitting errors to the individual probabilities of those age-specific survival data must be much less than the above. Hence, function (38) is the relation among the components of the age-specific survival vector we are looking for.
Practical identification algorithms Using the curve-fitting method of (39), agespecific survival data like the direct solution of (31) can be smoothed easily no matter how dispersive they are, and therefore the parameter estimation problem of the age-process can be solved. This, however, requires two steps: first, solve (31), and second, solve (39). Now substituting p/(at) in (40) for q in (1), we obtain a new model structure of the age-process which has only seven independent parameters. Substituting the q(=) for t/in (31), we can therefore rewrite the parameter estimation problem of the age-process as LefL scare .u U
Time-invariant and time-varying parameter estimation results of the age-process Assuming the age-specific survivals to be timeinvariant in the two periods, 1953-1964 and 1964-1975, respectively, and following the above method, we then obtain two sets of optimal parameter estimates of the age-process, q,(~) and r/z(~), which we call the time-invariant results of identification of the age-process. In order to obtain the dominant law of the agespecific survivals with respect to time, we noted the following facts: first, the total mortality of the Chinese population has, basically, been decreasing exponentially over the 25 years, with the exception of a few years of unusual mortality--the unusual part of which has been considered as an external disturbance; second, in a steadily developing society and under the condition that the senectitude states are basically stationary during the time-region under discussion a reduction in the total mortality must, to some extent, correspond to some reduction Rig.t
I
,...
r
scale
0 80
0 99
i, I,
I
"%~
Fitting 098 • o,
Fi%tea
t
060 ~-
o,~ e,,,. 0
2,0
40
60
(42)
instead of (31). It can be seen that the dimension of the nonlinear optimization problem is reduced from I01 to 7, and its computation-time is evidently reduced too. Using the DFP-algorithm without the analytic gradient, we can also do it on a DJS-130 computer. Since relation (38) represents the internal relation among the components of the age-specific survival vector q, it can be regarded as a further identification of the usual discrete model structure, and this is just the distinctive point of this paper in treating the identification problem. And that is why we call this kind of parameter estimation a parameter identification.
I
m
•
q
+ 2n + t Ile,{k + n, 1/(~t)}lt
80 Age (years)
I00
Fla. 3. A real example of fitting (381.
{20
140
The parameter identification of a population model of China of the age-specific mortalities. So we can assume the dominant law of the age-specific mortalities with respect to time to be exponentially decreasing, too; and ~h (~) and r/2(¢) obtained above to be equal to the age-specific survivals at the moments t~ and t2 = tl + n, respectively. Thus the general expression of the survival function can be written as r/(¢, t) = 'l ~.(¢) - c(~) exp { - :t(¢)t } ,1(¢.',tl ) = '11 (¢') ,7(¢,
(43)
t2) = ~2(¢)
where q=(~), c(~) and ~(~) are functions to be defined. Furthermore, we assume ~/=(~.) to be dependent on q:(~) and find the following form of ~(¢, t)
,ff~, t)
=
1 - k; + kcr/:(~) - [1 - k; + keq2(¢) - )h(~.)] 1 (1 -
k:)[1
-
"~.,+,)/.
r/=(¢)]
-
(44)
where k{ and t~ are both parameters to be determined under a definite time-origin. Thus the time-varying parameter estimation problem of the age-process can be roughly formulated as: choose the parameters k~ and t] so that the precision of total population forecast for each year in the past 25 years is optimal according to the survival parameters determined by (44). Thus we obtain the time-varying results of identification of the age-process which are expressed by the survival function (44) with respect to ~ and t, by means of which we can evaluate the age-specific survivals at any moment, not only from 1953 to 1978 but also after 1978. The former are called timevarying identification results and the latter the I /
parameter forecasting (or extrapolating) results corresponding to the former. When new statistical data come, we can use them to correct the earlier results, to obtain new results, and to make them more precise.
Validation of the identification results of the ageprocess Validity of the identification results can be verified only by studying the forecast precision of the past behaviours of the process, i.e. taking the population states in 1953 and 1964, respectively, as known initial distributions; the birth population every year as the inputs; and the parameters obtained above, time-invariant and time-varying, as those of the process model. We have calculated the population states in 1975, the total population and the total death population for each year in the 25 years, and compared them with real statistical data in Figs 4-6. Validation of the results with respect to both the time-invariant and the time-varying parameters obtained is basically satisfactory. The forecast precision of the total population for each year over the 25 years is within about 1 ~ . The validity and the extrapolation properties of the time-varying results are better than those of the time-invariant ones. So the forecast results of the parameters are valid too. IDENTIFICATION OF THE BIRTH-PROCESS
As mentioned in Appendix C, the identification of the birth-process here depends on the determinability of parameter estimates, or rather, it depends on whether the regression-solution of (33) or (34) is unique or not. Hence, if the time-varying parameter vector b(k) is considered as the direct object to be estimated, (33) is obviously not determinable. Even
J95~ -!975
:o
:
/~
v.:-
I
1964-1975 ~
°
k/
i
:o /),o ,i ,'o
'h//"l,"l
',1 'h
J964- 1975~
("~
(oI
1953-1975
~.ge (yeors ]
e,, -:0
423
/'~
I ' ~ I' " # I. I f % I/
(b)
FIG. 4. The predicting precision of the states in 1975: (a)by time-invariant results, (b)by time-varying results.
424
YUE-XIN ZHU and BAI-WU WaN
eee
*,eo
} Trine-varying result
...... I Time lrlvor~ont • • • • • A~ r e s u l t -
+I* Q
3
:
s
Q)
-S
t ,
\
io • ~ 3 ~ ~0
~58
,
• •~ • ~A73
,
__.&,
-,5 -- -
1953-1978
*
FIG. 5. The predicting precision of total population.
L 40
__-~____ } 1 9 5 3 - 1 9 7 8 • • • •
IA
r,'\--\
oI [ -'°b
',
',
..,.
..,-.m:" .o.
I ~.i¢ A" I *e~-'k I ~e 63T i ee r3o%~,,le78 ..... ........ " ' -
-e°/~-
. . . . ~ Time-varying • • • • • JI r e s u l t ...... Time-invoriant • • A• • result
F I G . 6. T h e p r e d i c t i n g
precision
of the total death
population.
the birth-process is mainly determined by the rank of X (see Appendix C). Obviously, if the population state is stationary (Kwakernaak, 1977), then, no matter how numerous the birth-statistical data are, the birthprocess is still not determinable. Although a theoretically stationary population does not really exist, a population process close to a stationary one, is generally not determinable either, because it usually requires finding the inverse of an illconditioned matrix. Thus, it follows that this kind of data-structure of the Chinese population system does not have universality in identifying the birthprocess, although it can be improved, to some extent, by means of a special algorithm. Fortunately, population states of China were not stationary in the 25 years and the rank of X, at least, is more than 10 which can be seen from its population states. In the Chinese population model, b(k) changes with respect to time mainly due to fl(k), while the fertility-pattern vector 2(k) is usually rather stable. So we use the following smooth function
b(~, t, v) = k/fll (t, v)21 (~, v) if b is assumed to be time-invariant, (34), and therefore the birth-process, is not determinable, either. This is because (rankX)m~=l<~-0+
1
where l = 25 as in (34); in other words, the birthstatistical data are too limited. However, as seen in identifying the age-process, the components of h are not independent of each other; a lot of real statistical results concerning the age-specific fertility-rates (for example, Ravenholt and Chao, 1974), show that there is some relation among them. In fact, nowadays some smooth functions are often used to express the dependent relation among the age-specific fertility-rates (Song and Yu, 1981). The most common one is the X2-probability density function
1 En/2- 1 e - ~ , 2 , 2,:2F(n/2) Px"(~) =
0
~< 0
e>0 (45)
where ~ = ~ - 0 + I; and n is a positive integer. This is used to describe the fertility-pattern function (Song and Yu, 1981), or rather, the relation among the components of the fertility-pattern vector 2, which is the same as b except for a scalar coefficient k:~ [see (8)]. On the other hand, it can be proved that if the relation among the components of b can be described by a smooth function with r independent parameters, and so long as rank X >/ r, (34) is then determinable. In the other words, with respect to our particular data structure, the determinability of
(46)
to identify the birth-process, in which k: is taken as a constant, equal to 0.4835 (Song and Yu, 1981); v' = [% vt...v9] is the independent parameter vector with dimension 10; while t - 8.005 fll(t,v) = Vo - vlt + v2tan - t - 0.001
exp{-O.5(t-v31V'~,t~[1953,1978], ~V,~l)
(47)
describes the time-varying properties of b and is determined by studying the change of the total birth-rate in the 25 years 21(¢, v) =
-
0 =, 1
'"
exp
-0.5
"V6
- 0 + 1
V8
(48) which describes the age-distribution of fertility-rate, or rather, the fertility-pattern, and is just a generalization of (45). In fact, it is the same as (45) provided that, in (48) v6=% =vo=
1 and
vT=n/'2- 1
except for the factor 1/[2"'2F(n/2)] which has been considered in (47). Thus, discretizing (48) and (47), according to (8) and (46) we can write the fertilitypattern vector 2' = [2o )-0~~... 2~] as
2 = 2t (v)/~/,.~ (a, v) where k't(v)= [).lt0, v)
,;.l(0+ 1, v)...21(~,v)]
(49)
The parameter identification of a population model of China
425
and the average fertility-rate of women as fl(k + j) = fl1(k + j,v) T. }q(a.v), a=O
j = O , 1..... 1 - 1.
(50) FIG. 7. The validation of the identification result of the birthpro~ss.
Substituting (49) and (50) into (33), we can write the linear regression equation (33) as the following nonlinear regression equation ~'
=
vector, obtained in this paper, as the model parameters, and by using the average fertility-rates of women, /3, as the parameter governed by the national family-planning policy. From Fig. 8, it follows that, in order to carry out the national population plan of China in which the total population is specified as under 1200 million by the end of this century, the average fertility-rates of women must be controlled under 1.70, and then the maximum of the total population will be about 1310 million and will take place in 2026 A.D.
DX).I (v)
where D = k/diag [fl~(k + 1 - 1, v) fl,(k + l -
2,v)... fll(k,v)].
So the parameter estimation problem of the birthprocess can be then converted into the nonlinear optimization problem ...,o
7 IIDX & (v) -
1
~'ll&.,
CONCLUSIONS
by which the birth-process parameter estimates, i.e. the fertility-pattern vector and the average fertilityrate of women, can be obtained. These results can be validated by the precision of calculations of the birth population in the 25 years, as shown in Fig. 7. Except for a few points, their precision is within 5 ~o-
APPLICATIONS
The parameter estimates obtained in this paper, along with the model, have been used to study prediction and control of the future trends of Chinese population in order to provide background material for the national population and familyplanning policy decisions. For example, Fig. 8 presents a few predictions for the total population, which were acquired by using the time-varying agespecific survival function and the fertility-pattern
From the above we conclude that the series of methods for identification seems to be effective in solving the parameter estimation problem of a population model under the conditions of very limited statistical data. This series includes implementation of a discrete model, decomposition and transformation of the bilinear model, a curve-fitting technique and nonlinear optimization, and, most important, the reduction of its dimensionality. The applicability condition for the curve-fitting technique is that there exists a functional relation among the parameters under estimation. The time-varying result makes the model more effective. This series of methods may be useful for population model identification in developing countries. Acknowledgements--The authors are indebted to Professor Bao~ sheng Hu, Huan-chen Wang, and Zhao-yong You for the discussions and their valuable suggestions.
t
B=2.20
g ,4 -._L3d_ _ _ . / ~ - e ~ - _ _ _
,o I.--
II
2026 6-
J980
2000
I
2020
I
2040
1
20¢~
I
~POQO
A.D. ( y e o r )
FIG. 8. A few predictions of total population of China by means of the identified parameters.
426
YUE-XIN ZHU and BAI-WU WAN
REFERENCES AIj, A. and A. Haurie (1980). Description and control of continuous-time population process. I E E E Trans. Aut. Control, AC-25, 3. Boyarsky, A. Ya and P. P. Shusherin (1955). Demographic Statistics. Gosstatizdat, Moscow (in Russian). Eykhoff, P. (1974). System Identification. John Wiley, London. Foerster, H. V. (1959). Some remarks on changing populations. In F. Stohman (Ed.), The Kinetics of Cellular Proliferation. Grune & Stratton, pp. 382-407. Isermann, R. (Ed.) (1979). Tutorials on system identification, 5th I F A C Symposium on Identification and System Parameter Estimation, Darmstadt. Keyfitz, N. (1968). Introduction to the Mathematics of Population. Addison-Wesley. Kwakernaak, H. (1977). Application of control theory to population policy. In A. Beusouson and J. Lions (Eds), New Trends in Systems Analysis. Springer, pp. 359-378. Leslie, P. H. (1945). The use of matrices in certain population mathematics. Biometrica, 33, 183. Ravenholt, R. T. and John Chao (1974). World fertility trends. Population Report, Family Planning Programs, Series J, No. 2. Song, Jian, and Jingyuan Y u (1981). On stability theory of population systems and critical fertility rates of women. Acta Automatica Sinica, 7, 1 (in Chinese). Wang. Huanchen (1980). Population state and population dynamics, Journal o f Xi'an Jiaotong University, 14, I (in Chinese). Zhu, Yue-xin (1981). Research on the parameter identification of a population model of China. Master's degree thesis, Xi'an Jiaotong University, China (in Chinese).
n-I
Hi,,(k)w(k + n - i - 1) = H~(k}w,(k).
(55)
i=O
Substituting (54) and (55) into (53), we then obtain (12a}. 2. Noting that yo(k + j ) = b'(k + j)xlk + j )
= b'(k +j){Hj(k)x(k) + H~(k)yoj(k) + H~(k)w~(k)}, j = 1, 2. . . . . n and b'(k +j)H~(k)yoi(k) = O, j = 1, 2. . . . . n < 0 a n d from (15) and (22a, b), we have H~(k)wl(k ) = H~a(k)wx(k),
(56)
j = 1, 2 . . . . . n.
Hence yo(k + j ) = b'(k + j)Hj(k)x(k) + b'(k + j ) H ~ (k)w,(k), j=l,
(57)
2,...,n.
Substituting it into (14), we have y6,(k) = [b'(k
I
I
I
I
'
1)H,_ t(k)x(k)l...~ ' b'tk)Ho(k)x(k)lOl.., i0]
+ n -
I
I
+ [h'lk + n - I)/-P'~ (k)w,(k)[...l
t
b'lk + 1)H~t(k)w,,(k) ~)0i...i 0]' I
I
I
j
•
I
+n-l),...i
+
w',(k) [ / ~
t'(k)b(k
+
n
-
Io]
1 )),' '" II
H~i(k)~(k + 1)i01. ..lO] ' A P P E N D I X A. A DERIVATION OF (12) A N D (13)
x'(k)B',,(k) + w',(k)B~'(k)
I. Letting u(k) = qo,,(k)yo(k) + w(k)
(51)
the transpose of which is just (13a)' 3. From (18), we have
from (9) we have
2,
x(k + 1) = H(k)x(k) + u(k)
sin(k)
x(k + 2) = H(k + l)H(k)x(k) + H(k + l)n(k) + n(k + 1)
=
,,,,k+2,.
__ ,,,,,,k.
L;(k+ n)
Li'x(k + n)J Li'H.(k)J
l'Ht(k)yo,(k)
tl-I
x(k + h) f Hh(k)x(k) + ~ Ho,(k)u(k + h - i -
x(k)+
=
Fl'H~(k)w,(k)]
1),
i=O
h = 0. I ..... n
(52)
where Hh(k), Hi.h(k) are given in (20) and
i'H~(k)yo,,(k).] Li'H•(k)w.(k)J Now, an arbitrary component of the second vector on the righthand side of (58) can be written as
-1
y [.] = o.
l'H~(k)yoh(k ) = q_~(k + h - l)yo(k + h - 1) + ...
i:0
+ q-I(k)qo(k + 1)...qh-2(k + h - l)yo(k)
Letting h = n and combining (51) and (52), gives
= l'H~l(k)yo,(k),
x(k + n) = H.(k)x(k) n-I
where H~l(k) is shown in (27b). Putting (59) and (56) into the second and third vector on the right-hand side of (58), respectively, we then obtain
+ ~ Hi.,,(k)qoa(k + n -- i -- 1)yo(k + n - i - 1) i=O ¢1-1
+ ~ H~.,(k)w(k+n-i-
1)'
(53)
i=O
s,(k) = F~(k)x(k) + G,(k)yo,(k) + L~(k)w,(k)
H~.,(k)qo,,(k + n - i - 1)yo(k + n - i - 1) =
]
Fd(k)
1
[ 0.p I I /I i=~'tr//(k + n - i +j)yo(k -~ n - i - 1)1~0 i . . . 1,~7 0 /j, _ i
i + t ..... m
Ld(k+n-l)J
and therefore + H~,.(k)qo.(k + n - i - l)yo(k + n - i - I ) = H~,(k)yo.(k).
F
(54)
L/~'(k+n-llx(k+n-1
{I -,lo,,(k)]yo(k) {! - r/o,(k + 1)}yo(k + 1)
LIt -
/=o
Noting (15)and (22a), we have
(60)
where Fs(k), G,(k) and L,(k) are (26a, b, c), respectively. Similarly, we have
From (2), (20), (3), and then (14) and (21), we have
0 ..... i - - 1
(59)
h = 1, 2 . . . . . n
l
,J (6l)
rlootk ÷ n - 1)}yo(k + n - 1)
where the vector~(k + j ) , j = 0, 1..... n - 1, is in (27a). The first vector on the right-hand side of (6l) represents the total death
The p a r a m e t e r identification of a population model of China population in state-transition process; and the second one the death number of babies born in the very year. Putting (12a), (27c) and (14) into (61), we have
d.(k) = / { (k + L'~'(k + n - 1)H._dk
+
because, at ~ ffi 0 u(t - ~)~.o = u(t) = p(O,t).
Hence the function
)]
x(~,t) =
po(~,t)d~,
~¢[-1,m-
x(a, k) = xo(k),
(67)
From (63)-(67), the function defined as ~¢[-1, m-1]
tp(~,t)=x(~+l,t+l)/x(~,t)
Substituting (591 and (56) into the above, we then obtain d.lk) = F~(k)x(k) + Galk)yo,,(k) + L,~(k)w,,(k)
(62)
where F.(k), Gd(k) and L,~(k) are in (26b), (26d) and (26f), respectively. Substituting (12a), (60) and (62) into (17), and noting (25), we have
/ = I r'(k) I x(k) ÷ l~'(k)lY°"(#)
Ld.(k) J L~',(k)_l
L~Ak)J
l
w.(k)l LL,(k) J |L,(k)
ffi Co(k)x(k) + D.(k)Yo,,(k) + Ro(k)w.(k)
which is just (12b). Equation (13h) is obvious. APPENDIX B. A PROOF OF THE PROPOSITION We assume it to be reasonable that continuous mathematical models are used to describe population process~. This is true, i n particular, when the total number of population individuals is large e n o u g h - - that is the way of Chinese population. In continuous population models (Focrster, 1959; Song and Yu, 1981), the state variable function is a continuous distribution function, with respect to
p(~,t)>o, ~ C0, m] which represents the population number per unit age at age ~ and at moment t; and the input variable function is a continuous distribution function, with respect to t
u(t)>O,
t~(-~,
+oo),
which represents the birth-population number per unit time at moment t, and, when ~ ffi 0
u(t) = p(O, t).
~
is just the age-specific survival function ~(~) at moment t we are looking for. This is because, except for x(~,t) = 0, the function tp(~, t) is obviously differentiable with respect to ~; and, for V~ ¢ [ - 1, m - 1 ], the function ,p(~, t), in the sense ofthe discrete model is just the survival at age ~ and at moment k; and, in particular, ~p(~,t)k=,=q.(k),
a = - l , O, .... m - 1 .
tlk
APPENDIX C. ON DETERMINABILITY OF PARAMETER ESTIMATES By means of the least-square method, the identification of both the age-process and the birth-process is transformed into a minimization problem. Obviously, especially under the condition of those very limited observed data, it is essentially important that whether the solution of the minimization problem, i.e. the parameters under estimation, be unique or not. If not, the process will not be identifiable anymore by those data. Conversely, if the solution is unique, it will be the most probable solution of the parameters determined by the observed data. Although some other arbitrariness still exists, we call this process determinable under the observed data, and regard this solution as the final result. 1. Definition. For certain observed data available, a process is said to be determinable under the observed data if its model parameters can be uniquely determined by those data. 2. Determinability o / t h e age-process. Determinability of the age-process under the data structure here depends on whether the optimal solution of (31) with respect to ~ is unique or not. However, we have Theorem 1. The optimal solution of the nonlinear optimization problem (31) with respect to ~ is unique. Proof. It is enough to prove that, without the statistical data i, a and ~,(k), this solution of (31) would still be unique. In fact, letting ~J.(~) 8q
we can then analYtically find the unique optimal solution qo as a-I
7° = ~,+1(k + n)/~o(k + n - a - 2 ) / r ]
which is just yo(k) in (3b) when above equation, we have
:o
qJ,
j=-I
t+ 1
yo(t) ffi
(68)
~o ffi :~o(k + n ) / ~ o ( k + n + 1)
Under this condition, the aggregate number of the birthpopulation born in [t,t + 1) can be written as
yo(t) =
(66)
a = O, 1. . . . . m
x ( - I, k) = yo(k).
°
+
(65)
and, from (63)
~ ' ( k + , - 1)n:-'Ik)yo.,- ~(k)J
z,.(k + n) |s.(k)
1]
is differentiable with respect to ~ and is the aggregate population number in [~, ~ + 1) at t. When ~ ffi a and t ffi k, then
~.'(k + 1)Hl(k)yodk)
I
427
a--0,1 ..... n-2
u(~)d~ a-1
t ffi k. L e t t i n g
~
--- t - ~ in the
~°ffi~.+1(k4n)l~.-.+,(k)/
H j-e--n+
~, 1
a=n-l,n-2~
u(t - Od~.
(63)
-1
We now define an auxiliary function
~
~E I-l,0)
(64)
which is obviously continuous with respect to ~ ¢ [ - 1 , m ]
.... m - l .
The detail of the proof is omitted here. 3. Determinability o f the birth-process. If the parameter vector h is taken as a direct object of the parameter estimation, the birth-process described in (33) will obviously be not determinable; but we have Theorem 2. The birth-process described in (34) is determinable if and only if rank X ffi ~ - 0 + I.
428
YUE-XIN Z r t c a n d B A b W U W A N
Proof. If
birth-process is rank X ~ ( - 0 +
1
rankX<~-0+
1
rank X ~> r.
there must be
Proof. Discretizing
rankX=~-0+l there must be either
v) b(O +
I, v), .... b((, v)].
Since b(~, v) is continuous and differentiable in expressed as
v-space, it can
be
(69)
b(v)=b(vo)+{~}'Av+O(Av)
(71)
I
rank [X II~'] = rank X or
b(~, v), we have
h'(v) = [b(0,
and (34) must have infinitely many solutions of b and is not determinable. Conversely, if
(70)
I
wherev : Vo + Av; 0(Av) is higher order infinitesimal vector than Av. Neglecting 0(Av) in (71), and substituting it into (34), we have
AAv
rank [X ll~ ] > rank X In the former case, (34) is consistent and has a unique solution of b. In the latter ease, (34) is contradictory and its regression solution b = (X'X)-
= g
(72)
where • f~b'(vo))
A= A ~ ;
'
(73)
tX'i,
exists uniquely because of (69). /..emma 1. The birth-process is not determinable if I < ( - 0 +1. /,emma 2. The birth-process of a stationary population under this data structure is not determinable. These two lernmas both are obvious because
g = ~ -
Xb(vo)
with dimension I × r and I × 1, respectively. For the vector Vo of arbitrary choice, the matrix A and the vector g can both be regarded as known, and the linear regression problem (34) with respect to b has been transformed into (72), with respect to Av. From Theorem 2, the neoessary and sufficient condition of the determinability of (72) is rank X = r. From the theorem of rank of a product matrix, this theorem is then proved. Remarks. 1. From the above proved, it follows that
( r a n k X)max = l < ~ - 0 + l
rank
and rankX = I respectively. 4. If b(~, v) is the relation among components of the timeinvariant vector b where v is the independent parameter vector with dimension r; and b(~, v) is continuous and differentiable in vspace; then we have:
Theorem 3.
A necessary condition of determinabifity of the
= r
(74)
is another necessary condition of determinability of the birthprocess. Because the condition (74) is often satisfied easily here, the condition (70) also has some suff~ientness, or rather, the determinability of the birth-process is mainly determined by (70). 2. When the relation among components of the time-varying parameter vector b(k) is expressed as (46), the conclusion about the determinabifity of the birth-process is the same as (70), and its proof is not presented here.