- Email: [email protected]
A Method for Modelling with Missing Observations
Copvright © IF.-\C Id entificat ion and S,'stelll Par;lllleter Estilllation 198:,. York . L't\.. I (I~F)
A METHOD FOR MODELLING WITH MISSING OBSERV ATIONS
,
Zhang Zhong-jun* and Wang Zhiz-hong** *DI'P(II'I III f//I of AlllolIWlioll. Shllllghlli J illo TOllg [ ' /lii'l'nilv. SIII/IIghll i. Thl' Pl'()pll' 's RI'IJ/liJ/i1' of' Chillll **Depll rl 1111'11 I of PrecisiulI 111.111'111111' 111. Shllllghlli Jillll Tllllg [·II;"l'nily. S//{/IIg/lIIi. Thl' Pl'IJpll''s Rl'plliJ/ic ofChill1l
Abstract. This paper proposes a new method for e~tabli~hin R discrete models with missing observations. The authors firstly develop a formula for calculating the unbiased estimator of the variance of noise serie~ {e,} ; then the model parameters are eAtimated by Monto-Carlo Atatistic method. This proposed method is suitable not only for establishing AR mOdels in time series analysis, but also for e~tablishing AHY~ models. For input-output descret ~ systems with missing observationR in input and output data. this method is also appplicable for modelling. .,:'
Keywords. l'arameter estimation; Missing observations; Hierarchical systems; Optimisation; ~lonte-Carlo methods.
1. INTRODUCTION (2.2)
J
To identify many practical systems, the required observed data may be incooplete. e.g. a part of sampled data had been missed due to one or the other following causes: 1. Observed data were not completely taken; 2. The sensor for recording data failed to operate properly; 3. Some data are unreliable due to intensive interference; 4. One sensor is used by time-sharing in measuring and record in" different time serie" data sequence". Quintana (1982) investigated the problem to estimate system states with unreliable observati.ons for electrical power systems. Tugnai t (1983) did reseach work on state estimation for systems with interrupted observationA. But their workA a~ mainly based on Rubstituting the conditional expectations at the missing points as the misAing data. Up to present still insufficient work has been ~evoted to expound the parameter estiMation problem for Rystems with missing data.McGiffin (1980.1981) proposed a method to eqtimate parameters ot· AR models with missing observations. The authors in this paper present a new parameter ~stimation method with missin6 data, which is applicable not only for establishing All models. but also for establishing AH!·:A models. ,'or input-output deAcrete systems with missing observations in input an~ output data. this method is also applicable for modelling.
the estimated parameters of the model or the estimator i\ =(b~ ..... b, .-a, ..... _a.)' is obtained from the following condition: mint J) = J
9=
(2.4) '
y,
y,
Y..,
y,.,
y,
,Vi
x
y", y
u ...... ,., U't"OI! Y.: Yi Ul-' - - - - - - - - - -r- - - - - - - - - U'."_"HI ••• U, .., J .Y;"t l
,
.
.v ... .
Yi+"I
Y,.".,
I: I
U ....H1
Y...
y-.~
(2.5) Howevpr. eq.(2.4) is not directly applicable for the case where some observation data are missing. To tackle this problem. let us develop the followin~ theorpm. On bases of thiq theoreM. we propose 3. new method. (Theorem) ::luppose y" .,Vi, •••• 'Yi, in {Y.) sequence are missing. with i"i ...... 1. as the serial numbers of the mi s sin"j data. then the unbiased estimator of the \'ariance of (e.) is
Conoider the following s ystem model
(2.6)
0.-'= S / (N-n-m-p) (2.1)
where
3 = min J
I e=e·'YI.=y,: .... ,y,,=y,; (2.7)
(proof) Define E = [e.., .p~, ••••• e~.1T , the n the reA1dual equation is
+ a,z"+ ••• + a.z·· , •• + b~z(n and m are assumed known)
B(Z")
XT)(f'xTy Uti
NOISE VARIANCE
A(Z") =
(2.3)
where
n. THE UNBIASED ESTHIATOR. Or'
where
I e. e
As le,} is an uncorre la ted ze ro mean noise sequence. 9 can be rea
1
= b, z"+
E
= Y - )(.A
and eq.(2.2) can be
tu,) and {Y.} are input ann out. put sequenc~s respectively. "hile le,) is an uncorrelhted zero mean noise ~eq uence. The number of samplec points is N+n. To estimate parameters of the model froOl ohserved data ~equcnces often requires to define a quadratic performance criterion dS follows:
(?R) ~ritten
as
J ~ ET],; (2.9) Obviously. S Cfln not he f'xpressf'o hv an explicit solution any longer. Becau~e all of the mis~ing data are included in X and y. the residual vf'ctor E can be E'Xprp.3 :,,~j by a nonlinear function
14~l7
14~H
Zhang Zhong·jun and Wang Zhiz·hong (2 .1 0 )
=j{
I::
T
whf'rp
~x ~ ansi~ B
f'. )
at t~ 2 vicinity of ~. into a Tay3p.rics. afti'r irooning a ll but the linpar
lo~
tt")rng, we havp.
r. ,;f ( ~· )
M(~)I 4>=4>' ' 4
:1.
T
whprp thp incrcnlf'ot
f ind th" ;n i ni:nwn
TO
,Yi,1
T
0
(2.11)
At:; = (..:18 T ,A Yi, t~ "{h f ••• ,AYi,1 T f .J '
p.qu~tion
ing itera.tin g
E(l(·')= Attc } . ~.IJ('+ ElM'
(2.12)
where (k) o r (k+1) r ep re sent t he iteration numbprn; thE' matrlx A=nf/rt~T has thn following f o rm
..
A
-a~ .... 1_'1e..... ~b,: '3a"
l l.... " 'c)b",
Thus
S
= min
J
.-
lim min ( F,' ... 'TE,,·.. ,)
2
'1 e ••"
-ae .... :lJc .....
-'3a, :'1Yi,
?Y"
.,.= 0
( 2.14)
Substituting e'1.(;>.12) into eq.(2.14). we have 2 Af.I'l T ( A ~J ,,4>"'U() + E (le)} = 0 '~""=- CA ,<) TA"')"'
1. ().
A,·IT E'"
( ;>.15)
Now E'··'ITE lll4") = (A(.).A(P"·~ ")+ r;lll'f (A'·) . A4' ..(tc J +
ing step. 'Cl ~'/ ' ", y .... ,=y,., _~_J_ ~ (' 4) nalculate v . ( i 1JY'rc, y .... retun to step (?). ' '" ,
El.']
a,r' A"" £··... E'·')
= l-A."'( A"lT A""'" A.""T £"'+E""J T l-A"'(A."IT A
(3.1)
Por calculating e'1.(3.1). let us neiluce the expres"ion for -aJ,onhyr" and TJ'"'hy.""'. For simplir:1ty. thp supersc~ipt" (k) are omitted in subsequent lines. J
(Y_X9)T (Y-Xq)
(;> .13)
In each iteratton the necessary condition for the existance of minimizing is gran (I::'.... TE'.. ·').,•
w,c:)••.••
seaching n irections t and take w/,". uJ,~·) as step factor". which equal the relaxation factor ,,'~devided by Recond partial d"rivatives ?'J"hy""l (i=i, ••••• i, ). From the above discussion WE' suggp st thp. follOwing hierarchical optil1li7.atio!'l alg orithm; (1 ) Gues'l initial value" for missing da ta y; ( i=i, ••••• i,). say y,=(y,.,+y, .. )/2 (when the nata Re'1upnCe has close co rrelatio n ,
yTy _ yTX(XTXr'XTy
(>.2)
Taking partial derivative of J w.r.t. y,. we have HT
~=
2 Y_2"y'X9_2 y TH 9+26T XTn 9 'by, "y, '/) y. ~y,
IiY,
(3.3)
Noting the position of ~ in X and Y in eq.(?5) whf'n n
U.4)
l 0 ••• 01; 0 ... 0: 0 ••• 0)· X9 (U;.~•••• Ui_, ; Yi-..... Y,_I ) • e
l1Y\9
,y
E/IOT(I._ A'·'(A,"'TA"')"'A'·JT)E'"' ::K
trace(I .. - A(IC'(AVIJ'T A(allf' AllIl ' ] r-: tll ) EliOT
hen~e
th~
(;>.16)
B(z")U,-(1 - A(Z"))y, ~
('xPE'ct at i on of 3 i a
lim
E7~race
U!'
{
lim {N
17"' -
I I 'f'dA
"'T
9 X -
'ilY,
(2.17)
A
9 =
eT
0
si (N-n-l1I-p)
I
ui
IY'
Y.. -,""
In the las t ~pction. we sep. that if the minimum of eq.(2.7) is founo. thp.n th" :mbiaspd efltimator .. 'ca n be ohtn inpd from e'1.(2.17). HOllever. a " lobal optimizatio n algorithm to seak the conver5ence point ~. in eq.(?7) is a time-consuming task. In t his Rp.ction. «e propose a t,.(o If'vel hip.rarc hir:al opti miza tion a l gorithm for calculating 1;hp. solution o ~· p.o.(2.7). Eac h missing datuo is taken as a "pconn l<>vel subs.ystpm; minimize J with respect to (w.r.t.) parameter vector e aB the fi rst level. 1'ill;.3.1 shows tne structure of this algorithm.
I
:
I
•
IY' MI.I
(2. Hl)
Ill. THE HIERARGHICAL O?'lU!IZATION AWORITHM lOR ESTIr-:t.TlN(; THE NOISE VARIANCE
.
I U i .".,
Lp.. the unbiased esti!!lato r of cr' is
&-' '"
(3.6)
u, .•• ,
IU,,.,,.,,,,
(1.- A'"'(A''''A,IIfr'A'·'T)·
(N-n-m-p)
(3.5)
[ o... 0 i Y•••••• y, •. J . 9
E (.3) = lil1l E { I-: ""'T £ ""'}
.-
yt
e
Y,
( o... 0 ; y.t.... y,!. ). El
(3.7)
Jubstituting eqs.(3.4)-(3.7) into eq.(5.3). we have
~~ '3
= 2y. i
2yt -
2f{0 •.• O:y •••••• Y... )
'f
f"
- (0 •.. O;y,...... y ... )}9
2e~ - 2
r o ... 0 i ef..... et.. J . e
..
2(et + a.e~ + ••• + ~here
a.e!..)
7.A (z)e~
(3 .6)
et = y. - yf
(3.9)
In the s ame way . we g~t • Zt- •• • +a " ........ z ••• ' .) E",• ~ ={ <.~ (1 Ta.
'3y~
2 (A. ... ,., Z"ooi •• + •• .
+a. z·
In ord p. r to minimize J. wP. may take minus partial derivatives of J w.r.t. each missing data as the
) e~
N
(5.10)
" .U'.. l
J,.
fi g.) .' Thp hierarchical optimazation alg ori th"1 for minimizing J
:\ Methud fur Modelling with l\lissing OiJsen"ations ferent narametf'r es timat o r s . (4) C~lc ' llatinp, thp uv " mgc value s of par:.mpter estimator" ohtained in stop (5), thp.y ~ill be taken 303 thp pardmeter estimators of the monel.
+
8o,,(y, .., +8.,
"
Y;:tlH
+ •.. +a y; ....-b. u i 9t
+_.- .. .
=2(1+a~+ •.• +a:)
.
.!!...
D.11)
.
dimilarly
, -b",u, ... ~J
" + •.. +a • ) = { 2( 1+a,
(3.12)
2(9. ..\ .. + •.• +8.:)
'aY,l
discussion is o~ly for the situation data. For input sequence {u,) with missing data, we have similar expressions: Th~
In thp case ·...hpre bo th \Y.) and {u.) have Jli. Rslng ohservations, we JGa y morllfy the aho ve rcont o-:a r lo mpthorl : Firstly, for in put by arl~in G white noi s p scquence ~ith pr0per i ntesity to ~ " ... ,u~ till J=(N-n-m-p)~', thp.n for Olltput by arloiing another noi~e spquence to Yi , ' " ' ' Yi, till J=( N-n-m l iT-'.
abov~
~ith out~ut mis~in0
(j=j"j., ••• ,j,)
(3.13)
V. The genp-raliza tion AR we kno~, paran,gter esti;n"tion by usine GLS (Genp.mli7.ed Lgast dquares) is more accurate. Lpt the mo~el identified by ilL:'; haR the following :'orrn 1 ('i.' ) A( 7. .') Yt= B ( z .') u t + C(:>:") ct where A(z") and B(7.·') havp the same forms af> in
eq.(?.1); C(z·' )=1+c,7.·'+ ... +c,z·'. Put y,=C(z ") y, ann li.. =C(z·')u,, then eq.(5.') woulrl h~ rewritten as In a similar ~ay w~ can gpt similar hierarchical optimization ftlgorithm. The only rlifference is ~J
1Y ·
;<1 '3J
-=L--=- ' ~ 11 y, j.; 11 Yj 'Hi
1iJ ?J = - +c-+
, "aY",
raY~
By usinl~ the ah"vp- I)xpressions and thp. hierarchical optimi7.ation algorithm, ~. can be found readily. In this case ~'= J/(N-n-m-p-q).
=
1J~
-ay.
(' +c, Z + ••• +c, z, )
= 2G(z)·A(7.)e! I V. 'rHE ~10N'rO-CAlU.O ImTl!OV r'OR PARAr:'ETi.1
ESTIMATION ',HTH
~n3SING
DATA
With the hierarchical optimization algorithm. ~e can get not only the unbiased p.stimator of ~'. but also thp. ~. ~hich satisfies eo.(2.7). It must be emphasized that the vector e· ohtained ahove is not yp-t the model parameter estimator. AS ~how ahovp-. the expectation of J is (N-n-m-p)~l, while the expectation of J of the original model without missing data iA (N-n-m)~·. The statistical characterA of these t~o models are different. AO e· obtained above cann't be considered as the unbiased estimator of e. In fact, for a stochastic process it is impossible to determine the exact values of missing data, "hat we can do most is to establish a model from observation sequences with thp- same statistical properties without missing data. In order to maintain the same second order moment of noise sequence, we suggest the MontoCarlo metho~ as follOWS: Add zero mp-an white noise with proper intensity to Yi,'''''Yi, obtained above. Att/>=.p", J hal' the minimum value. Thus by adding noil'e with any intensity to Yi,' ••• ,Y." J ~ill bp. certainly increased. Increasing J from (N-n-m-p)~'up to (N-n-m)~' USing searching method with appropriate white noise repeatedly, we get several slightly different paramp.tp.r e ,nimators. The averaee valup.s of them may b~ t ~ ken as the parameter ~Rtimator of the model. Jwnmari7.ing thp above discussions. '.re state our 8ugges tpd methon for estimating model parameters with missinb observa tions as follo~s: (1) Calculatp ~. and y~ •••• ,y~ hy u~inF, hi~rar chical o~timi7.ation al30rithm. (2) Add 11 z? ro mean p:;euno-mndoJ~ ~hit~• .n0i3e sequencp, of ~hich the variance equals k ~ • to Y.~ ••. , ,ft, . By seaching, I
(5.3)
Evioently, if w~ Rubstitute C(Z)A(Z) ann C(z)B(z) for A(Z) and B("z) in sec.IlI. a Similar method can h~ oht'lir.erl. The Rimula tion rpsul ts aN' also sa tisfactory. To p.stiQate parameters of an AR model with missing observations the procedure is just the same aR in spctions II-IV except B(z~)=O in eq.(2.1). However. in the case of ARI'iA mortels. the calculation will be more complicated. Assume the following AHMA monel
~
~ - D(z~) e,
where
(5.4)
C(z~)=1+c.z~+
D(z-')=1+d,
--
+C"'Z
z~+
As we know, an A~~ model can be equivalent to a an AR model ",i th infinite terms, ?8 expressp.d by the following equation A( 7..,
whpre
)r, = et
(5.5)
A(z")=1+a, z"+a.z··+ ... + ...
::)uppose <#1= (c ....... ,Cl ,-0", ... I-d, 'Yi,,· .. ,:r;,)" ,
thp-n the above-mentioned method can be applicable for pammetpr eiltimation of ARI';}' m0dels with consideration of the following aspects: ( 1) El in an ARMA mod"l has to be c" lcula ted h.Y using nonlinear leaRt s'1uares, such ail ;1arquardt algorithm (1963). (2) Aftpr minimzing the performance criterion J of an AR mod~l with missing 1ata. take the Y~"" •.• ,y~ of this AR monel a ~ the i~itial v"lu~A of these miSSing data in the ~ip.rarchic a l optimization al,o;ortthm of the AIWEA morlpl. so ail to l"Pouce the numh"r 0" c",llin.'O the ,w"routine of nonlinear leas t "~\.L·H·~S for rerlllc in,:; thp comp 'lta tion time. D) '3J/tay, and 1i'J/~y.. are 1iffprent fro", tho1lC! in A~ model3. ~or an infinite AR modpl, ( <; .7)
Zhang Zhong-jun and Wang Zhiz-hong
I"O() Com;,ini:lF; wit " "q.(').6).
~
h'lve
','?
VI. CONCLU,lIONS
.Q~z.l. ~(Z") 'Hi - 2 C(z) C~z") Yi
_
How€":~r
~lJ/1J'y.l~
(,).8)
a.2
~
= 1 +
J'"
In thi" p~per a new llnthoi for efltablishin~ Aft and inp ,lt-out:)c.t "YRtpm monels is oror:oo",,,rl. r~(' C:1mnl..l.t r r ,~il:".ulation rpsillts are BRtiRfar.-
(').9)
AXi'.A
J
AR Box (1976) llpntion"a. '.tlpn thf' model satiRfie" invprtibility ~onctit~on (Le, . pvpry roo ts of C(w)=O ,,"'9 l oo:·,tf'd Ollt of trf' 'lnit circle) ~, '>() A.nri
thprl 3
HellOP
t
Um ['l1J /'tIY: l..,.ce t....,..o
-(1+J~'
REFERENCES
a'l )
Urn (g,e"!l· f'·'Y(1_e"g,)]
I; ,,,-2&.j=
Box, G.E.P., and G.M.Jnnkins (1176). Time serips analysis forf'castin" and control. rev.ed., !:lan ~'rtmcisco, Calif .• /!olden-Day. ~iarQuardt (131):5). An algoritlur. for l,.,ailt-sqLl."res estimation of nonlinear pa!'ameterR. J. of SIAN. D. ~31-441. --McGiffin, P.B .• anrl D.N.i'.Hurthy (19 nO,1981). i'arar'lp.te r e3timation for auto-reGressive sys tems wi th mis3 in;; obse rva tions -- pa rt I ann part 11. Int. J. of Ryst. sCi., 11. 10211034. and 12, (,5 7-663. -~uint?na. ·{.H. (19R2). Had data n,.,tection anrl identification techniques usine estimation orthogonal ml'tholis. lEEr; trans., pas-l01, 55'5(,-536~. --'fugnait, J.K. (19,9 3). On identification and adaptiv<' f'"ti'na tion for sy"tpms with interruntect ohsprvations. Automatica, 19. 6173. -
t.,. ...
.Fltl
< c.. (3. small p03itive r.umber) (5.11) ThuR we may sele'cot a sufficiently 19.r" e numher 1 to get. " satiAf""to"'y app roximate value of'liJ/'iJY:. a, •...• iL in nq.(').'1) can hp ohtained by 1ividing il(.,) '::!y C("'),
'I. ;3H:ULATIOt; RESULTJ
l EXclmp1e 1 J Usp a Zf'ro mean Garssian white noise sequence 9." (e.) . Among 100 sampled data. abnut 1/3 is miflsinp,. U"in~ non1inear leaAt square". to efltimate the narnmeters of the follow in,; AHnA model
=( 1 +0. 5 £' ) "'t
( 1- 1 • ~, £' +0 . 7 Z·, ):It
,
~ABr,II 1 Rpsul+' o"'irlulat:ion of ARI'lit. ~10de1 +
"
'"
er
d,
d,
c,
True values
1.0
-1.5
0.7
0.5
Without missing data
0.9fl3
-1.513
0.719
0.494
Wi th periodic missing data
0.982
-1.539
With missing data ___.
0.886
-1.542
O.;;J
r -- "random L
I
0.:4~
0.570
[V:xa!Jlnle 2) USf' two un~orrel"ted 7.f'ro mean ;}au~sian white noisD sequence w~th variances 1 and (j' respectiv,.,ly as {u.} and { Po. ) • i.. monp, lOO sampled ~ata, ahout 1/4 is mifl"ine. Using JL!:l to estimatp. th" parametern of the following model (1 -1.5£' +1. 4z" - O.fl5z-').v. = (z-' -z" +0. 5z")u. +e.
I( 1-0. 5z-' ... 0. 3z-')
TABLS 2 Results of !:limulation of Input-Output Model
True values Without missing data Cl! +' '+-I..,
'"
co I, e!l
h~
o
III ~CJ
....a
I .... 0
s..Ql.., ....
P-oO
..,s:: ~
I.... 0
~
L~
0",
r-.
....ID E
!
.... Cl.., ~
P-o 0
..,cE 0
~
a,
a.
b,
b.
bs
1.4
-0.85
1.0
-1.0
0.5
0.1026 -1.499 1.402 -0.858 0.994 -0.999 0.497 0.1061 -1.507 1.416 -0.861 0.995 -1.0230.535 0.1024 -1.501 1.406 -0.859 0.989 -1.003 0.506
{u..) &: {Y.}
0.1031 -1.505 1.410 -0.860 0.994 -1.012 0.519
( Ut)
0.1005 -1.501 1.408 -0.861 0.983 -0.999 0.505
{ Y.)
0.1031 -1.502 1.409 -0.860 0.993 -1.004 0.511
I u.)
&:
/Y.}
0.1058 -1.506 1.410 -0.858 0.992 -1.021 0.524 1.0
-1.5
'.4
-0.85
1.0
-1.0
0.5
1.035
-1. 366 1.350 -0.781 0.918 -0.941 O. 364 ~ -1.430 1.412 -0.809 0.865 -1.04') 0.708
/Y.)
1.021
-1.414 1.390 -0.797 0.934 -1.029 0.554
{uJ &: (Yt)
1.047
-1.403 1.386 -0.7850.791 -0.965 0.439
{Ut}
0.976
-1.399 1.384 -0.807 0.042 -0.981 0.411
{Y.}
1.012
-1.384 1.362 -0.778 0.889 -0.943 0.390
(u.) &: {Y.}
1.055
-1.405 1.404 -0.780 0.861 -').955 0.512
Without missing data I .... 0
a, -1.5
(Y. )
True values
Cl! +'
u0.1
{u.}
a0
(u.)
1.016
obsprvatiollf1 de lose in-
,"orn or 1",,,,, . Evirlpnt1y, t.he l~ss missin< Doint" in the o~servef. "equencp.. the hetter fo,.. lToorlt"lling.
(C;.10)
8,j') = U"lli t...,.o<> j"'l1"1
Ti:-;~iIlg
f OrnJat~ on
ga'>O j>l
"Urn ~
t or-y .. HO,",f!ver,
A METHOD FOR MODELLING WITH MISSING OBSERV ATIONS
,
Zhang Zhong-jun* and Wang Zhiz-hong** *DI'P(II'I III f//I of AlllolIWlioll. Shllllghlli J illo TOllg [ ' /lii'l'nilv. SIII/IIghll i. Thl' Pl'()pll' 's RI'IJ/liJ/i1' of' Chillll **Depll rl 1111'11 I of PrecisiulI 111.111'111111' 111. Shllllghlli Jillll Tllllg [·II;"l'nily. S//{/IIg/lIIi. Thl' Pl'IJpll''s Rl'plliJ/ic ofChill1l
Abstract. This paper proposes a new method for e~tabli~hin R discrete models with missing observations. The authors firstly develop a formula for calculating the unbiased estimator of the variance of noise serie~ {e,} ; then the model parameters are eAtimated by Monto-Carlo Atatistic method. This proposed method is suitable not only for establishing AR mOdels in time series analysis, but also for e~tablishing AHY~ models. For input-output descret ~ systems with missing observationR in input and output data. this method is also appplicable for modelling. .,:'
Keywords. l'arameter estimation; Missing observations; Hierarchical systems; Optimisation; ~lonte-Carlo methods.
1. INTRODUCTION (2.2)
J
To identify many practical systems, the required observed data may be incooplete. e.g. a part of sampled data had been missed due to one or the other following causes: 1. Observed data were not completely taken; 2. The sensor for recording data failed to operate properly; 3. Some data are unreliable due to intensive interference; 4. One sensor is used by time-sharing in measuring and record in" different time serie" data sequence". Quintana (1982) investigated the problem to estimate system states with unreliable observati.ons for electrical power systems. Tugnai t (1983) did reseach work on state estimation for systems with interrupted observationA. But their workA a~ mainly based on Rubstituting the conditional expectations at the missing points as the misAing data. Up to present still insufficient work has been ~evoted to expound the parameter estiMation problem for Rystems with missing data.McGiffin (1980.1981) proposed a method to eqtimate parameters ot· AR models with missing observations. The authors in this paper present a new parameter ~stimation method with missin6 data, which is applicable not only for establishing All models. but also for establishing AH!·:A models. ,'or input-output deAcrete systems with missing observations in input an~ output data. this method is also applicable for modelling.
the estimated parameters of the model or the estimator i\ =(b~ ..... b, .-a, ..... _a.)' is obtained from the following condition: mint J) = J
9=
(2.4) '
y,
y,
Y..,
y,.,
y,
,Vi
x
y", y
u ...... ,., U't"OI! Y.: Yi Ul-' - - - - - - - - - -r- - - - - - - - - U'."_"HI ••• U, .., J .Y;"t l
,
.
.v ... .
Yi+"I
Y,.".,
I: I
U ....H1
Y...
y-.~
(2.5) Howevpr. eq.(2.4) is not directly applicable for the case where some observation data are missing. To tackle this problem. let us develop the followin~ theorpm. On bases of thiq theoreM. we propose 3. new method. (Theorem) ::luppose y" .,Vi, •••• 'Yi, in {Y.) sequence are missing. with i"i ...... 1. as the serial numbers of the mi s sin"j data. then the unbiased estimator of the \'ariance of (e.) is
Conoider the following s ystem model
(2.6)
0.-'= S / (N-n-m-p) (2.1)
where
3 = min J
I e=e·'YI.=y,: .... ,y,,=y,; (2.7)
(proof) Define E = [e.., .p~, ••••• e~.1T , the n the reA1dual equation is
+ a,z"+ ••• + a.z·· , •• + b~z(n and m are assumed known)
B(Z")
XT)(f'xTy Uti
NOISE VARIANCE
A(Z") =
(2.3)
where
n. THE UNBIASED ESTHIATOR. Or'
where
I e. e
As le,} is an uncorre la ted ze ro mean noise sequence. 9 can be rea
1
= b, z"+
E
= Y - )(.A
and eq.(2.2) can be
tu,) and {Y.} are input ann out. put sequenc~s respectively. "hile le,) is an uncorrelhted zero mean noise ~eq uence. The number of samplec points is N+n. To estimate parameters of the model froOl ohserved data ~equcnces often requires to define a quadratic performance criterion dS follows:
(?R) ~ritten
as
J ~ ET],; (2.9) Obviously. S Cfln not he f'xpressf'o hv an explicit solution any longer. Becau~e all of the mis~ing data are included in X and y. the residual vf'ctor E can be E'Xprp.3 :,,~j by a nonlinear function
14~l7
14~H
Zhang Zhong·jun and Wang Zhiz·hong (2 .1 0 )
=j{
I::
T
whf'rp
~x ~ ansi~ B
f'. )
at t~ 2 vicinity of ~. into a Tay3p.rics. afti'r irooning a ll but the linpar
lo~
tt")rng, we havp.
r. ,;f ( ~· )
M(~)I 4>=4>' ' 4
:1.
T
whprp thp incrcnlf'ot
f ind th" ;n i ni:nwn
TO
,Yi,1
T
0
(2.11)
At:; = (..:18 T ,A Yi, t~ "{h f ••• ,AYi,1 T f .J '
p.qu~tion
ing itera.tin g
E(l(·')= Attc } . ~.IJ('+ ElM'
(2.12)
where (k) o r (k+1) r ep re sent t he iteration numbprn; thE' matrlx A=nf/rt~T has thn following f o rm
..
A
-a~ .... 1_'1e..... ~b,: '3a"
l l.... " 'c)b",
Thus
S
= min
J
.-
lim min ( F,' ... 'TE,,·.. ,)
2
'1 e ••"
-ae .... :lJc .....
-'3a, :'1Yi,
?Y"
.,.= 0
( 2.14)
Substituting e'1.(;>.12) into eq.(2.14). we have 2 Af.I'l T ( A ~J ,,4>"'U() + E (le)} = 0 '~""=- CA ,<) TA"')"'
1. ().
A,·IT E'"
( ;>.15)
Now E'··'ITE lll4") = (A(.).A(P"·~ ")+ r;lll'f (A'·) . A4' ..(tc J +
ing step. 'Cl ~'/ ' ", y .... ,=y,., _~_J_ ~ (' 4) nalculate v . ( i 1JY'rc, y .... retun to step (?). ' '" ,
El.']
a,r' A"" £··... E'·')
= l-A."'( A"lT A""'" A.""T £"'+E""J T l-A"'(A."IT A
(3.1)
Por calculating e'1.(3.1). let us neiluce the expres"ion for -aJ,onhyr" and TJ'"'hy.""'. For simplir:1ty. thp supersc~ipt" (k) are omitted in subsequent lines. J
(Y_X9)T (Y-Xq)
(;> .13)
In each iteratton the necessary condition for the existance of minimizing is gran (I::'.... TE'.. ·').,•
w,c:)••.••
seaching n irections t and take w/,". uJ,~·) as step factor". which equal the relaxation factor ,,'~devided by Recond partial d"rivatives ?'J"hy""l (i=i, ••••• i, ). From the above discussion WE' suggp st thp. follOwing hierarchical optil1li7.atio!'l alg orithm; (1 ) Gues'l initial value" for missing da ta y; ( i=i, ••••• i,). say y,=(y,.,+y, .. )/2 (when the nata Re'1upnCe has close co rrelatio n ,
yTy _ yTX(XTXr'XTy
(>.2)
Taking partial derivative of J w.r.t. y,. we have HT
~=
2 Y_2"y'X9_2 y TH 9+26T XTn 9 'by, "y, '/) y. ~y,
IiY,
(3.3)
Noting the position of ~ in X and Y in eq.(?5) whf'n n
U.4)
l 0 ••• 01; 0 ... 0: 0 ••• 0)· X9 (U;.~•••• Ui_, ; Yi-..... Y,_I ) • e
l1Y\9
,y
E/IOT(I._ A'·'(A,"'TA"')"'A'·JT)E'"' ::K
trace(I .. - A(IC'(AVIJ'T A(allf' AllIl ' ] r-: tll ) EliOT
hen~e
th~
(;>.16)
B(z")U,-(1 - A(Z"))y, ~
('xPE'ct at i on of 3 i a
lim
E7~race
U!'
{
lim {N
17"' -
I I 'f'dA
"'T
9 X -
'ilY,
(2.17)
A
9 =
eT
0
si (N-n-l1I-p)
I
ui
IY'
Y.. -,""
In the las t ~pction. we sep. that if the minimum of eq.(2.7) is founo. thp.n th" :mbiaspd efltimator .. 'ca n be ohtn inpd from e'1.(2.17). HOllever. a " lobal optimizatio n algorithm to seak the conver5ence point ~. in eq.(?7) is a time-consuming task. In t his Rp.ction. «e propose a t,.(o If'vel hip.rarc hir:al opti miza tion a l gorithm for calculating 1;hp. solution o ~· p.o.(2.7). Eac h missing datuo is taken as a "pconn l<>vel subs.ystpm; minimize J with respect to (w.r.t.) parameter vector e aB the fi rst level. 1'ill;.3.1 shows tne structure of this algorithm.
I
:
I
•
IY' MI.I
(2. Hl)
Ill. THE HIERARGHICAL O?'lU!IZATION AWORITHM lOR ESTIr-:t.TlN(; THE NOISE VARIANCE
.
I U i .".,
Lp.. the unbiased esti!!lato r of cr' is
&-' '"
(3.6)
u, .•• ,
IU,,.,,.,,,,
(1.- A'"'(A''''A,IIfr'A'·'T)·
(N-n-m-p)
(3.5)
[ o... 0 i Y•••••• y, •. J . 9
E (.3) = lil1l E { I-: ""'T £ ""'}
.-
yt
e
Y,
( o... 0 ; y.t.... y,!. ). El
(3.7)
Jubstituting eqs.(3.4)-(3.7) into eq.(5.3). we have
~~ '3
= 2y. i
2yt -
2f{0 •.• O:y •••••• Y... )
'f
f"
- (0 •.. O;y,...... y ... )}9
2e~ - 2
r o ... 0 i ef..... et.. J . e
..
2(et + a.e~ + ••• + ~here
a.e!..)
7.A (z)e~
(3 .6)
et = y. - yf
(3.9)
In the s ame way . we g~t • Zt- •• • +a " ........ z ••• ' .) E",• ~ ={ <.~ (1 Ta.
'3y~
2 (A. ... ,., Z"ooi •• + •• .
+a. z·
In ord p. r to minimize J. wP. may take minus partial derivatives of J w.r.t. each missing data as the
) e~
N
(5.10)
" .U'.. l
J,.
fi g.) .' Thp hierarchical optimazation alg ori th"1 for minimizing J
:\ Methud fur Modelling with l\lissing OiJsen"ations ferent narametf'r es timat o r s . (4) C~lc ' llatinp, thp uv " mgc value s of par:.mpter estimator" ohtained in stop (5), thp.y ~ill be taken 303 thp pardmeter estimators of the monel.
+
8o,,(y, .., +8.,
"
Y;:tlH
+ •.. +a y; ....-b. u i 9t
+_.- .. .
=2(1+a~+ •.• +a:)
.
.!!...
D.11)
.
dimilarly
, -b",u, ... ~J
" + •.. +a • ) = { 2( 1+a,
(3.12)
2(9. ..\ .. + •.• +8.:)
'aY,l
discussion is o~ly for the situation data. For input sequence {u,) with missing data, we have similar expressions: Th~
In thp case ·...hpre bo th \Y.) and {u.) have Jli. Rslng ohservations, we JGa y morllfy the aho ve rcont o-:a r lo mpthorl : Firstly, for in put by arl~in G white noi s p scquence ~ith pr0per i ntesity to ~ " ... ,u~ till J=(N-n-m-p)~', thp.n for Olltput by arloiing another noi~e spquence to Yi , ' " ' ' Yi, till J=( N-n-m l iT-'.
abov~
~ith out~ut mis~in0
(j=j"j., ••• ,j,)
(3.13)
V. The genp-raliza tion AR we kno~, paran,gter esti;n"tion by usine GLS (Genp.mli7.ed Lgast dquares) is more accurate. Lpt the mo~el identified by ilL:'; haR the following :'orrn 1 ('i.' ) A( 7. .') Yt= B ( z .') u t + C(:>:") ct where A(z") and B(7.·') havp the same forms af> in
eq.(?.1); C(z·' )=1+c,7.·'+ ... +c,z·'. Put y,=C(z ") y, ann li.. =C(z·')u,, then eq.(5.') woulrl h~ rewritten as In a similar ~ay w~ can gpt similar hierarchical optimization ftlgorithm. The only rlifference is ~J
1Y ·
;<1 '3J
-=L--=- ' ~ 11 y, j.; 11 Yj 'Hi
1iJ ?J = - +c-+
, "aY",
raY~
By usinl~ the ah"vp- I)xpressions and thp. hierarchical optimi7.ation algorithm, ~. can be found readily. In this case ~'= J/(N-n-m-p-q).
=
1J~
-ay.
(' +c, Z + ••• +c, z, )
= 2G(z)·A(7.)e! I V. 'rHE ~10N'rO-CAlU.O ImTl!OV r'OR PARAr:'ETi.1
ESTIMATION ',HTH
~n3SING
DATA
With the hierarchical optimization algorithm. ~e can get not only the unbiased p.stimator of ~'. but also thp. ~. ~hich satisfies eo.(2.7). It must be emphasized that the vector e· ohtained ahove is not yp-t the model parameter estimator. AS ~how ahovp-. the expectation of J is (N-n-m-p)~l, while the expectation of J of the original model without missing data iA (N-n-m)~·. The statistical characterA of these t~o models are different. AO e· obtained above cann't be considered as the unbiased estimator of e. In fact, for a stochastic process it is impossible to determine the exact values of missing data, "hat we can do most is to establish a model from observation sequences with thp- same statistical properties without missing data. In order to maintain the same second order moment of noise sequence, we suggest the MontoCarlo metho~ as follOWS: Add zero mp-an white noise with proper intensity to Yi,'''''Yi, obtained above. Att/>=.p", J hal' the minimum value. Thus by adding noil'e with any intensity to Yi,' ••• ,Y." J ~ill bp. certainly increased. Increasing J from (N-n-m-p)~'up to (N-n-m)~' USing searching method with appropriate white noise repeatedly, we get several slightly different paramp.tp.r e ,nimators. The averaee valup.s of them may b~ t ~ ken as the parameter ~Rtimator of the model. Jwnmari7.ing thp above discussions. '.re state our 8ugges tpd methon for estimating model parameters with missinb observa tions as follo~s: (1) Calculatp ~. and y~ •••• ,y~ hy u~inF, hi~rar chical o~timi7.ation al30rithm. (2) Add 11 z? ro mean p:;euno-mndoJ~ ~hit~• .n0i3e sequencp, of ~hich the variance equals k ~ • to Y.~ ••. , ,ft, . By seaching, I
(5.3)
Evioently, if w~ Rubstitute C(Z)A(Z) ann C(z)B(z) for A(Z) and B("z) in sec.IlI. a Similar method can h~ oht'lir.erl. The Rimula tion rpsul ts aN' also sa tisfactory. To p.stiQate parameters of an AR model with missing observations the procedure is just the same aR in spctions II-IV except B(z~)=O in eq.(2.1). However. in the case of ARI'iA mortels. the calculation will be more complicated. Assume the following AHMA monel
~
~ - D(z~) e,
where
(5.4)
C(z~)=1+c.z~+
D(z-')=1+d,
--
+C"'Z
z~+
As we know, an A~~ model can be equivalent to a an AR model ",i th infinite terms, ?8 expressp.d by the following equation A( 7..,
whpre
)r, = et
(5.5)
A(z")=1+a, z"+a.z··+ ... + ...
::)uppose <#1= (c ....... ,Cl ,-0", ... I-d, 'Yi,,· .. ,:r;,)" ,
thp-n the above-mentioned method can be applicable for pammetpr eiltimation of ARI';}' m0dels with consideration of the following aspects: ( 1) El in an ARMA mod"l has to be c" lcula ted h.Y using nonlinear leaRt s'1uares, such ail ;1arquardt algorithm (1963). (2) Aftpr minimzing the performance criterion J of an AR mod~l with missing 1ata. take the Y~"" •.• ,y~ of this AR monel a ~ the i~itial v"lu~A of these miSSing data in the ~ip.rarchic a l optimization al,o;ortthm of the AIWEA morlpl. so ail to l"Pouce the numh"r 0" c",llin.'O the ,w"routine of nonlinear leas t "~\.L·H·~S for rerlllc in,:; thp comp 'lta tion time. D) '3J/tay, and 1i'J/~y.. are 1iffprent fro", tho1lC! in A~ model3. ~or an infinite AR modpl, ( <; .7)
Zhang Zhong-jun and Wang Zhiz-hong
I"O() Com;,ini:lF; wit " "q.(').6).
~
h'lve
','?
VI. CONCLU,lIONS
.Q~z.l. ~(Z") 'Hi - 2 C(z) C~z") Yi
_
How€":~r
~lJ/1J'y.l~
(,).8)
a.2
~
= 1 +
J'"
In thi" p~per a new llnthoi for efltablishin~ Aft and inp ,lt-out:)c.t "YRtpm monels is oror:oo",,,rl. r~(' C:1mnl..l.t r r ,~il:".ulation rpsillts are BRtiRfar.-
(').9)
AXi'.A
J
AR Box (1976) llpntion"a. '.tlpn thf' model satiRfie" invprtibility ~onctit~on (Le, . pvpry roo ts of C(w)=O ,,"'9 l oo:·,tf'd Ollt of trf' 'lnit circle) ~, '>() A.nri
thprl 3
HellOP
t
Um ['l1J /'tIY: l..,.ce t....,..o
-(1+J~'
REFERENCES
a'l )
Urn (g,e"!l· f'·'Y(1_e"g,)]
I; ,,,-2&.j=
Box, G.E.P., and G.M.Jnnkins (1176). Time serips analysis forf'castin" and control. rev.ed., !:lan ~'rtmcisco, Calif .• /!olden-Day. ~iarQuardt (131):5). An algoritlur. for l,.,ailt-sqLl."res estimation of nonlinear pa!'ameterR. J. of SIAN. D. ~31-441. --McGiffin, P.B .• anrl D.N.i'.Hurthy (19 nO,1981). i'arar'lp.te r e3timation for auto-reGressive sys tems wi th mis3 in;; obse rva tions -- pa rt I ann part 11. Int. J. of Ryst. sCi., 11. 10211034. and 12, (,5 7-663. -~uint?na. ·{.H. (19R2). Had data n,.,tection anrl identification techniques usine estimation orthogonal ml'tholis. lEEr; trans., pas-l01, 55'5(,-536~. --'fugnait, J.K. (19,9 3). On identification and adaptiv<' f'"ti'na tion for sy"tpms with interruntect ohsprvations. Automatica, 19. 6173. -
t.,. ...
.Fltl
< c.. (3. small p03itive r.umber) (5.11) ThuR we may sele'cot a sufficiently 19.r" e numher 1 to get. " satiAf""to"'y app roximate value of'liJ/'iJY:. a, •...• iL in nq.(').'1) can hp ohtained by 1ividing il(.,) '::!y C("'),
'I. ;3H:ULATIOt; RESULTJ
l EXclmp1e 1 J Usp a Zf'ro mean Garssian white noise sequence 9." (e.) . Among 100 sampled data. abnut 1/3 is miflsinp,. U"in~ non1inear leaAt square". to efltimate the narnmeters of the follow in,; AHnA model
=( 1 +0. 5 £' ) "'t
( 1- 1 • ~, £' +0 . 7 Z·, ):It
,
~ABr,II 1 Rpsul+' o"'irlulat:ion of ARI'lit. ~10de1 +
"
'"
er
d,
d,
c,
True values
1.0
-1.5
0.7
0.5
Without missing data
0.9fl3
-1.513
0.719
0.494
Wi th periodic missing data
0.982
-1.539
With missing data ___.
0.886
-1.542
O.;;J
r -- "random L
I
0.:4~
0.570
[V:xa!Jlnle 2) USf' two un~orrel"ted 7.f'ro mean ;}au~sian white noisD sequence w~th variances 1 and (j' respectiv,.,ly as {u.} and { Po. ) • i.. monp, lOO sampled ~ata, ahout 1/4 is mifl"ine. Using JL!:l to estimatp. th" parametern of the following model (1 -1.5£' +1. 4z" - O.fl5z-').v. = (z-' -z" +0. 5z")u. +e.
I( 1-0. 5z-' ... 0. 3z-')
TABLS 2 Results of !:limulation of Input-Output Model
True values Without missing data Cl! +' '+-I..,
'"
co I, e!l
h~
o
III ~CJ
....a
I .... 0
s..Ql.., ....
P-oO
..,s:: ~
I.... 0
~
L~
0",
r-.
....ID E
!
.... Cl.., ~
P-o 0
..,cE 0
~
a,
a.
b,
b.
bs
1.4
-0.85
1.0
-1.0
0.5
0.1026 -1.499 1.402 -0.858 0.994 -0.999 0.497 0.1061 -1.507 1.416 -0.861 0.995 -1.0230.535 0.1024 -1.501 1.406 -0.859 0.989 -1.003 0.506
{u..) &: {Y.}
0.1031 -1.505 1.410 -0.860 0.994 -1.012 0.519
( Ut)
0.1005 -1.501 1.408 -0.861 0.983 -0.999 0.505
{ Y.)
0.1031 -1.502 1.409 -0.860 0.993 -1.004 0.511
I u.)
&:
/Y.}
0.1058 -1.506 1.410 -0.858 0.992 -1.021 0.524 1.0
-1.5
'.4
-0.85
1.0
-1.0
0.5
1.035
-1. 366 1.350 -0.781 0.918 -0.941 O. 364 ~ -1.430 1.412 -0.809 0.865 -1.04') 0.708
/Y.)
1.021
-1.414 1.390 -0.797 0.934 -1.029 0.554
{uJ &: (Yt)
1.047
-1.403 1.386 -0.7850.791 -0.965 0.439
{Ut}
0.976
-1.399 1.384 -0.807 0.042 -0.981 0.411
{Y.}
1.012
-1.384 1.362 -0.778 0.889 -0.943 0.390
(u.) &: {Y.}
1.055
-1.405 1.404 -0.780 0.861 -').955 0.512
Without missing data I .... 0
a, -1.5
(Y. )
True values
Cl! +'
u0.1
{u.}
a0
(u.)
1.016
obsprvatiollf1 de lose in-
,"orn or 1",,,,, . Evirlpnt1y, t.he l~ss missin< Doint" in the o~servef. "equencp.. the hetter fo,.. lToorlt"lling.
(C;.10)
8,j') = U"lli t...,.o<> j"'l1"1
Ti:-;~iIlg
f OrnJat~ on
ga'>O j>l
"Urn ~
t or-y .. HO,",f!ver,