- Email: [email protected]
Complete modeling of multidimensional noise as a multivariate autoregressive process
Ann. nucl.Energy,Vol. 18,No. 12,pp. 697-704,1991 Printedin GreatBritain.All rightsreserved
0306-4549/91$3.00+0.00 Copyright© 1991PergamonPressplc
COMPLETE MODELING OF MULTIDIMENSIONAL NOISE AS A MULTIVARIATE AUTOREGRESSIVE PROCESS N.
MORISHIMA
D e p a r t m e n t of N u c l e a r E n g i n e e r i n g , K y o t o U n i v e r s i t y , Yoshida, Sakyo-ku, Kyoto 606, Japan
( I August
1991)
A b s t r a c t - - A new m e t h o d for m u l t i v a r i a t e a u t o r e g r e s s i v e modeling of a vector time series is described, which has an a d v a n t a g e in g u a r a n t e e i n g l e s s - b i a s e d e s t i m a t i o n in a l e a s t - s q u a r e s sense by e n s u r i n g o r t h o g o n a l i t y between the c o m p o n e n t s of a residual v e c t o r . An expression for a m u l t i v a r i a t e a u t o r e g r e s s i v e model with a z e r o - l a g coefficient m a t r i x is derived t h e o r e t i c a l l y and utilized for the m o d e l i n g based on a weighted least-squares procedure. Model p a r a m e t e r s such as regression m a t r i c e s and residual m a t r i x are d e t e r m i n e d recurslvely for each value of o r d e r by using the LWR a l g o r i t h m . It is shown t h a t the present modeling differs f r o m the o r d i n a r y one in e s t i m a t i n g o p e n - l o o p p r o p e r t i e s , t h o u g h these m o d e l i n g are equivalent for c l o s e d - l o o p e s t i m a t e s .
1.
INTRODUCTION
For a vector time series { x ( n ) } to be linear and s t a t i o n a r y , m u l t i v a r i a t e a u t o r e g r e s s i v e ( M A R ) modeling is a c c o m p l i s h e d by the idea t h a t { x ( t t ) } is g e n e r a t e d f r o m a vector w h l t e - n o i s e process { r ( n ) } with m u t u a l l y independent components. We m a y thus model i x ( n ) } as the o u t p u t of a l i n e a r filter with coefficient m a t r i c e s A ( m ) in response to { ~ ( n ) } , t h a t is M
x(n)
-
-~,~(m)x(n-m)
+ 7(n),
(1)
m=l
where }*! is model o r d e r and ~: =
(r (n)}
.
has a yah'lance m a t r i x (2)
Note t h a t the s u p e r s c r i p t T denotes a t r a n s p o s e o p e r a t i o n . All the model p a r a m e t e r s such as A ( m ) and F. a r e easily d e t e r m i n e d by the l e a s t - s q u a r e s fitting of Eq. ( 1 ) to m u l t i d i m e n s i o n a l noise d a t a , since well tried and highly efficient c o m p u t i n g a l g o r l s m s a r e now a v a i l a b l e . H o w e v e r , the m a i n p r o b l e m is t h a t the c o m p u t e d residual m a t r i x E is not a l w a y s g u a r a n t e e d to be significantly dlagonal, t h a t is, the desirable modeling to obtain unbiased 697
698
N. MORISHIMA
e s t i m a t e s in a l e a s t - s q u a r e s sense is not a l w a y s g u a r a n t e e d . One a p p r o a c h to this p r o b l e m is to select a s a m p l i n g period so s m a l l as to yield a d i a g o n a l E(Klelss, 1982 ; Morishima, 1986) . A n o t h e r is to i n t r o d u c e a z e r o - l a g regression matrix A(0) w i t h the e l e m e n t s d e t e r m i n e d by a p r i o r i p h y s i c a l information(Kanemoto, 1982). The p r e s e n t p a p e r describes a new m e t h o d of M A R m o d e l i n g w h i c h a l w a y s yield s p e c t r a l e s t i m a t e s to be l e s s - b l a s e d in a l e a s t - s q u a r e s sense. This is a s s u r e d by a r e s i d u a l m a t r i x E d i a g o n a l i z e d at e a c h value of ~/, t h a t is, by m u t u a l i n d e p e n d e n c e of r e s i d u a l - n o i s e c o m p o n e n t s . This m e a n s t h a t the m o d e l i n g is the c o m p l e t e d e c o m p o s i t i o n of a vector time series as a n M A R p r o c e s s , w h i c h is a c c o m p l i s h e d by the f o l l o w i n g i d e a . F i r s t the r e s i d u a l m a t r i x E which is r e a l a n d s y m m e t r i c is f a c t o r i z e d into the p r o d u c t of o r t h o g o n a l m a t r i x P as E=PPT. T h e n E q . ( 1 ) is m u l t i p l i e d by P-I f r o m the l e f t - h a n d side. E x p r e s s i n g the result as a r e s p o n s e of { x ( n ) } to P - i v ( n ) , we thus obtain M
z(n)
=
~,[E~..0-A(m)}z(n-m)
+ r(n),
(3)
where p--1 ,
A(.~) =
p_l~(m),
f o r m= 0 for In _> 1
(4)
and
r(n)
= P-t-r(n),
(5)
w h e r e E is a unit m a t r i x a n d 8~.0 the K r o n e c k e r ' s delta defined to be u n i t y if ra=0 a n d zero if m # 0 . Obviously a new r e s i d u a l vector r ( n ) satisfies the r e q u i r e m e n t of m u t u a l i n d e p e n d e n c e for the c o m p o n e n t s since
( r ( n ) r ( n ) T) = p - t ( r ( n ) r ( n ) T ) p - T = p-1Ep -T
6) Another
f e a t u r e c a n be seen f r o m
(r(n)Tr(n))
=
the r e s i d u a l
(r(n)TE-'r(n)),
sum
of s q u a r e s 7)
w h e r e E -1 is served as a m a t r i x of w e i g h t s . Hence the above new m o d e l i n g is based on the w e i g h t e d l e a s t - s q u a r e s f i t t i n g of the generalized M A R model in E q . ( 3 ) , while an o r d i n a r y M A R m o d e l i n g is the case of a unit m a t r i x of weights a n d A ( 0 ) - E . The p r e s e n t p a p e r is o r g a n i z e d as f o l l o w s . C h a p t e r 2 is devoted to a t h e o r e t i c a l d e r i v a t i o n of the M A R model in E q . ( 3 ) as a g e n e r a l r e p r e s e n t a t i o n of a vector t i m e series. C h a p t e r 3 describes a m e t h o d a n d p r o c e d u r e for the modeling. It is s h o w n t h a t the p a r a m e t e r s A(m)(11~=O,1,2...M) a n d E c a n be d e t e r m i n e d r e c u r s i v e l y f o r e a c h value of /~1 by using c o n v e n t i o n a l a l g o r i s m s such as the LWR o n e . C h a p t e r 4 is concerned with s p e c t r a l p r o p e r t i e s of the e s t i m a t e s in t e r m s of c l o s e d - a n d o p e n - l o o p c h a r a c t e r i s t i c s . In p a r t i c u l a r ,
Multivariate autoregressivemodeling
699
it is shown that the magnitudes of residual variances can be chosen freely by adjusting the weight m a t r i x without a n y change in spectral estimates. Fi nal l y, c h ap ter 5 gives some concluding r e m a r k s .
2.
MAR REPRESENTATION
For a multivariable system to be linear and s t a t i o n a r y , consider a vector time series { x ( n ) } as a set of observations made at equidistant time intervals n h ( n = O , 1 , 2 . . . . ). In order to derive an MAR representation of { x ( n ) } , we first define a p o w e r - s p e c t r a l - d e n s i t y ( P S D ) m a t r i x P ( f ) for { x ( n ) } :
p(f)
(8)
zx ( X ( f ) X ( f ) * ) ,
where f is a f r e q u e n c y , X ( f ) is a Fourier component vector for { x ( n ) } , and the superscript t denotes the operation of H erm i t i an transposition. Since P ( f ) is H e r m i t i a n , it can be factorized into the diagonal m a t r i x E ( f ) of real positive eigenvalues and the o r t h o g o n a l m a t r i x U ( f ) of the eigenvectors:
P(f) We m a y rewrite Eq.
= U(f)E(f)U(f)' (9)
P(f)
(9)
as follows.
= U(f)E(f)~/2E(fN)-~/2E(fbl)E(fN)-~/zE(f)x/2U(f) ' = T ( f , f N ) - I E ( f N ) [ T ( f , f N ) - x ] *,
(10)
where fN is the Nyquist frequency equal to 1/2h and
T(f,fN)
= [U(f)E(f)I/2E(fN)-X/2} -I
(11)
Since T ( f , fN)-I is causal and asym pt ot i cal l y stable, it satisfies the condition of invertibility for a moving average process ( H a y k i n and Kesler, 1979). Hence we can express T ( f , f N ) in the form of an infinite Fourier series :
T(f,fN)
= ~-~,A(m)z -'h,
(12)
m=O
where z=exp(i27zf) given by
and A(m) is the Fourier coefficient m a t r i x of the expansion
A(m)
= f_tNT ( f fN
fN)z'hdf '
(m=0
1 2 ...). '
(13)
'
Hence we have expressed P ( f ) in terms of A(ra) and E ( f N ) . This leads us to represent { x ( n ) } as an MAR process in the form of Eq. ( 3 ) where the regression m a t r i x A(la) is defined by Eq. ( 1 3 ) and the residual-noise
700
N. MORISHIMA
vector r(Tt) has noise processes to be white and m u t u a l l y independent such that
(r(n)r(n')r>
=
Eti.,~,,
(14)
with the diagonal m a t r i x of variance fN = ./-/fE(fN)df
=
E(fN)/h.
(15)
Notice t h a t , taking the Fourier t r a n s f o r m of Eq. ( 3 ) and inserting the result on X ( f ) into Eq. ( 8 ) , we then obtain P ( f ) in Eq. ( 1 0 ) . Obviously the MAR expression thus derived is characterized by the z e r o - l a g regression m a t r i x A ( 0 ) and the diagonal m a t r i x P. of residual noise. This is in co n tr as t to an o r d i n a r y MAR model since the latter can be derived from Eq. ( 3 ) as follows. By s ubt r a c t i ng x ( a ) from the both sides of Eq. (3 and multiplying the result by A ( 0 ) from the left, we then obtain Eq. (1 where A(m)
= A(0)-la(m)
Y(n)
= A(0)-~r(n).
(m=l,2 ....
),
(16
and (17
Consequently the o r d i n a r y MAR modeling does not distinguish a z e r o - l a g correlation f r o m the others. This is one of the reasons why we often encounter an undesirable significant c or r e l a t i o n of residual noises in spite of a considerably short sampling period.
3.
MAR MODELING
The Y u l e - W a l k e r equation for the present MAR modeling can be derived as follows. Using Eq. ( 1 3 ) for A(m) and also the covariance function m a t r i x R ( l ) given by
R(l)
we m a y calculate Kesler, 1 9 7 9 )
the
A(m)R(l-,~
=
f_/~ S~ P ( f ) z -th d f
convolution sum
(l=0
of A(la)
'
-el and
'
+2
'
R(l)
= (i27rh)-'E(fN)f{T(f,fN)-'}*z'-ldz.
...), as ( H a y k i n
(18) and
(19)
•=0
Since {T(f,fN)-~} * is anal yt i c on and inside the unit circle in the z plane, we can p e r f o r m the integral in E q . ( 1 9 ) and arrive at the Yule-Walker
Multivariateautoregressivemodeling
701
equation
~, A ( m ) R ( l - m ) • =o
-
{ E{A(0)T]-~ 0
(l=O), (l>1).
(20)
Once the values of R ( l ) are given, all the model parameters A(m) and E can be determined uniquely by using Eq. ( 2 0 ) . To show this explicitly, we m a y rewrite Eq. ( 2 0 ) for an order M in the form n An(O)R(I)An(O)' - -~,An(m)An(O)-1"An(O)R(t-/n)An(O)' m=l
(l=1,2
....
(21)
~/),
and n
EM-
An(O)R(O)An(O)'
+ ~,An(m)An(O)-t.AM(O)R(rn)'AM(O) '
(22)
Note that these set of equations are similar in form to the o r d i n a r y Yule-Walker equations, provided that the values of AM(O)R(I)AM(O)* are given to solve An(re)AM(O) -1 and En. That is, once an adequate initial guess for An(0) and a set of R ( l ) ( l = O , 1 , 2 . . . M ) are given, we can calculate all the values of A M ( m ) ( m = I , 2 , . . M ) and EM. Such a calculation will be repeated for various values of AM(0) until a diagonal form of EM is obtained. This means that an iterative least-squares fitting is necessary for unbiased estimation o£ the model p a r a m e t e r s . However, using the following procedures, we can determine all the model parameters at once, not iteratively. ( 1 ) Define the o r d i n a r y Yule-Walker equations which are derived from Eqs. ( 2 1 ) and ( 2 2 ) : M
R(l)
= -)--~,An(m)R(/-m)
( / = 1 , 2 . . . . M),
(23)
m=l
and M
~M = R ( 0 )
+ ~--~,AM(m)R(-ra),
(24)
lI~=1
where
AM(m) = AM( O)-tAM(m) ,
(25)
EM = AM( 0 )-IEM( An( 0 )-1 ) ,
(26)
and
( 2 ) Solve Eqs. ( 2 3 ) and (24) to determine the values A M ( m ) ( m = l , 2 . . . . ~q) and EM. ( 3 ) Factorize En into eigenvalues m a t r i x D and eigenvectors m a t r i x U:
EM = UDU'. (4) obtain
Notice
the
correspondence
of
(27) between
Eqs.
(26)
and
(27),
we thus
702
N. MORISHIMA
[:M = D and also,
f r o m Eq.
and
AM(O) = U ' ,
(28)
(25),
aM(m)
ffi A M ( 0 ) ~ M ( m ) ,
(m=I,2
....
H).
(293
There is nothing special for the p a r a m e t e r calculation, since an o r d i n a r y LWR a l g o r i t h m is utilized as it is. Only a care should be taken of the factorization of EM. T hat is, the m a t r i x U must be determined in correspondence to the structure of EM. For e x a m p l e , in the case of the correlated residuals where the cross correlation~ E M ( i , j ) and E M ( j , i ) are significant, U must have eigenvector components in the elements U ( i , i ) , U(i,j), U ( j , i ) and U ( j , j ) . Such a m a t r i x structure of U can be constructed easily and uniquely, for e x a m p l e , by the use of J a c o b i ' s method for o r th o g o n al t r a n s f o r m a t i o n of a real s y m m e t r i c m a t r i x .
4.
PROPERTIES
OF
SPECTRAL
ESTIMATES
4. I. Closed-Loop Properf~ The present modeling is equivalent to the o r d i n a r y one for spectral estimation of o v e r a l l ( c l o s e d - l o o p ) p r o p e r t y . This is ascertainable in terms of a PSD m a t r i x . The Fourier t r a n s f o r m of Eq. ( 3 ) yields
P(f)
~ (X(f)X(f)¢) M
M
- ( Yq,AM(m)
) -'r.M ( EAM( m)'
M
=
)
~=0
I~=0
M
(E+E]M(m)z-'h) -' ~M (E+Y]~M'(m)z'h) -' • =1
Obviously the last result in Eq. MAR modeling.
(30)
M=I
(30)
is just the PSD given by the o r d i n a r y
4 . 2 . Open-Loop Property Distinct difference with r e g a r d to the PSD is the power contributions from each residual noise. This can be seen f r o m Eq. ( 3 0 ) in terms of the a u t o - P S D s due to each diagonal of F.M for the present modeling, and the those f r o m each diagonal of ~M for the o r d i n a r y modeling. An i m p o r t a n t point is that the f o r m e r preserves the total power of a signal by the sum of power contributions, in contrast with the latter neglecting all the power contributions f r o m the o f f - d i a g o n a l s of ~M. The other o p e n - l o o p characteristics such as a direct response f r o m one signal to an o t he r and a noise source PSD are also different between the present and o r d i n a r y MAP~ modeling. This will be shown in the following. Taking
Multivariate autorcgressivemodeling
703
the Fourier t r a n s f o r m of Eq. ( 3 ) and splitting the regression m a t r i x diagonal one A d ( f ) and an o f f - d l a g o n a l one A o ( f ) as
-EA/~(m)z -'h -
Ad(f)
into a
(31)
+ Ao(f),
m=0
we thus obtain the i n p u t - o u t p u t relation for X ( f ) : X(f)
where G ( f )
G(f)X(f)
(32)
+ F(f),
is a f r e q u e n c y - r e s p o n s e function m a t r i x given by C(f)
and F ( f )
-
-
-Ad(f)-*Ao(f),
(3 3)
is a noise-source vector defined by the a u t o - P S D m a t r i x O(f)
-
(F(f)F(f)'}
-
Ad(f)-lEM(Ad(f)-') ' .
Consequently, G ( f ) and Q ( f ) the values of A~/(ra) and EM.
4.3. Residual-Noise
(34)
depend on the method of MAR modeling t hrough
Maonitudes
As shown in Eq. ( 6 ) , the new MAR modeling yields a unit m a t r i x of residual, while the corresponding one in Chapter 3 uses the ei genm at ri x D as residual m a t r i x E,~. This suggests t ha t the magnitudes of residual noise m a y be chosen freely without a n y change in spectral estimates. To prove this, we will factorize I:• in the f o r m , instead of Eq. ( 2 7 ) ,
(35)
EM ffi UV*/~WV'nU ' ,
where V and W are under the condition
the diagonal
matrices
with
arbitrarily-chosen
D = VIn~fV I n .
(36)
T hen , f r o m the correspondence between Eqs. new p a r a m e t e r s as F.M = W so that the other same way as Eq.
and
diagonals
(26)
and
(35),
7~M( O ) = V-~nU * ,
parameters A(m)(m=l,2 .... (29). Hence we obtain
/q)
we m a y define
(37)
can
be determined in the
H
E A ~ / ( m ) z -'h = V - ' h A d ( f ) m=0
+ V-~nAo(f),
(38)
704
N. MORISHIMA
Consequently the frequency-response function m a t r i x G ( f ) PSD m a t r i x Q ( f ) are given by G(f)
and the noise-source
= -(V-~/2Ad(f ) )-~(V-t/2Ao(f) ) -
-Ad(f)-~Ao(f),
(39)
and Q(f)
= ( V - t / 2 A d ( f ) )-tlc( ( V - ~ n A d ( f ) ) - ~ ) ' -
Ad(f)-~O(Ad(f)-~) '
(40)
Obviously ~ ( f ) and Q ( f ) are the results given by Eqs. ( 3 3 ) and ( 3 4 ) , respectively, and therefore independent of the factorizatlon of ~t~ into a diagonal m a t r i x E• and an orthogonal one A ~ ( 0 ) .
5.
CONCLUDING
REMARKS
For a linear s t a t i o n a r y vector time series, the MAR model has been derived theoretically in the generalized form with a zero-lag regression matrix and independent residual noises. It has been shown that the MAR modeling can be performed easily using the o r d i n a r y algorithm for calculating recurslvely a set of model parameters up to a specified order. The present modeling would be useful in estimating open-loop characteristics, notably when a residual-nolse correlation is significant. For example, the combined use of the present and o r d i n a r y modeling would be promising for checking the effect of residual-noise correlation on open-loop estimates. Such a spectral analysis is in progress now and the results will be reported somewhere in near future.
REFERENCES
H a y k i n , S. and S. Kesler (1979) Predictlon-Error Filtering and M a x i m u m - E n t r o p y Spectral Estimation, in N o t t l i n e a r ~ e t h o d s o f S p e c t r a l A a a l ! l s i s , edited by S. H a y k i n . S p r i n g e r - V e r l a g : Berlin, pp. 9-72. Kanemoto, S, Y. Andoh, F. Y a m a m o t o , K. Kitamoto and K . Nunome ( 1 9 8 2 ) Identification of Pressure Control System Dynamics in BWR Plant by Multivariate Autoregressive Modeling Technique. J. Nucl. Sci. Technol. 19, 5 8 - 6 8 . Kleiss, E. B. ( 1 9 8 2 ) Determination of BWR Incore Power Feedback Effects by Radial Coherence Measurements. P r o q . N u c l . E n e r g y . 9, 6 6 5 - 6 7 4 . Morishlma, N ( 1 9 8 6 ) On the Frequency Resolution of Conventional and Autoregressive Spectral Estimates of Stationary Reactor Noise. Ana. Nucl. En. 13, 125-130. Wiggins, R. A. and E. A. Robinson (1965) Recursive Solution to the Multichannel Filtering Problem. J. Geophy. Res. , 70, 1885-1891 .
0306-4549/91$3.00+0.00 Copyright© 1991PergamonPressplc
COMPLETE MODELING OF MULTIDIMENSIONAL NOISE AS A MULTIVARIATE AUTOREGRESSIVE PROCESS N.
MORISHIMA
D e p a r t m e n t of N u c l e a r E n g i n e e r i n g , K y o t o U n i v e r s i t y , Yoshida, Sakyo-ku, Kyoto 606, Japan
( I August
1991)
A b s t r a c t - - A new m e t h o d for m u l t i v a r i a t e a u t o r e g r e s s i v e modeling of a vector time series is described, which has an a d v a n t a g e in g u a r a n t e e i n g l e s s - b i a s e d e s t i m a t i o n in a l e a s t - s q u a r e s sense by e n s u r i n g o r t h o g o n a l i t y between the c o m p o n e n t s of a residual v e c t o r . An expression for a m u l t i v a r i a t e a u t o r e g r e s s i v e model with a z e r o - l a g coefficient m a t r i x is derived t h e o r e t i c a l l y and utilized for the m o d e l i n g based on a weighted least-squares procedure. Model p a r a m e t e r s such as regression m a t r i c e s and residual m a t r i x are d e t e r m i n e d recurslvely for each value of o r d e r by using the LWR a l g o r i t h m . It is shown t h a t the present modeling differs f r o m the o r d i n a r y one in e s t i m a t i n g o p e n - l o o p p r o p e r t i e s , t h o u g h these m o d e l i n g are equivalent for c l o s e d - l o o p e s t i m a t e s .
1.
INTRODUCTION
For a vector time series { x ( n ) } to be linear and s t a t i o n a r y , m u l t i v a r i a t e a u t o r e g r e s s i v e ( M A R ) modeling is a c c o m p l i s h e d by the idea t h a t { x ( t t ) } is g e n e r a t e d f r o m a vector w h l t e - n o i s e process { r ( n ) } with m u t u a l l y independent components. We m a y thus model i x ( n ) } as the o u t p u t of a l i n e a r filter with coefficient m a t r i c e s A ( m ) in response to { ~ ( n ) } , t h a t is M
x(n)
-
-~,~(m)x(n-m)
+ 7(n),
(1)
m=l
where }*! is model o r d e r and ~: =
(r (n)}
has a yah'lance m a t r i x (2)
Note t h a t the s u p e r s c r i p t T denotes a t r a n s p o s e o p e r a t i o n . All the model p a r a m e t e r s such as A ( m ) and F. a r e easily d e t e r m i n e d by the l e a s t - s q u a r e s fitting of Eq. ( 1 ) to m u l t i d i m e n s i o n a l noise d a t a , since well tried and highly efficient c o m p u t i n g a l g o r l s m s a r e now a v a i l a b l e . H o w e v e r , the m a i n p r o b l e m is t h a t the c o m p u t e d residual m a t r i x E is not a l w a y s g u a r a n t e e d to be significantly dlagonal, t h a t is, the desirable modeling to obtain unbiased 697
698
N. MORISHIMA
e s t i m a t e s in a l e a s t - s q u a r e s sense is not a l w a y s g u a r a n t e e d . One a p p r o a c h to this p r o b l e m is to select a s a m p l i n g period so s m a l l as to yield a d i a g o n a l E(Klelss, 1982 ; Morishima, 1986) . A n o t h e r is to i n t r o d u c e a z e r o - l a g regression matrix A(0) w i t h the e l e m e n t s d e t e r m i n e d by a p r i o r i p h y s i c a l information(Kanemoto, 1982). The p r e s e n t p a p e r describes a new m e t h o d of M A R m o d e l i n g w h i c h a l w a y s yield s p e c t r a l e s t i m a t e s to be l e s s - b l a s e d in a l e a s t - s q u a r e s sense. This is a s s u r e d by a r e s i d u a l m a t r i x E d i a g o n a l i z e d at e a c h value of ~/, t h a t is, by m u t u a l i n d e p e n d e n c e of r e s i d u a l - n o i s e c o m p o n e n t s . This m e a n s t h a t the m o d e l i n g is the c o m p l e t e d e c o m p o s i t i o n of a vector time series as a n M A R p r o c e s s , w h i c h is a c c o m p l i s h e d by the f o l l o w i n g i d e a . F i r s t the r e s i d u a l m a t r i x E which is r e a l a n d s y m m e t r i c is f a c t o r i z e d into the p r o d u c t of o r t h o g o n a l m a t r i x P as E=PPT. T h e n E q . ( 1 ) is m u l t i p l i e d by P-I f r o m the l e f t - h a n d side. E x p r e s s i n g the result as a r e s p o n s e of { x ( n ) } to P - i v ( n ) , we thus obtain M
z(n)
=
~,[E~..0-A(m)}z(n-m)
+ r(n),
(3)
where p--1 ,
A(.~) =
p_l~(m),
f o r m= 0 for In _> 1
(4)
and
r(n)
= P-t-r(n),
(5)
w h e r e E is a unit m a t r i x a n d 8~.0 the K r o n e c k e r ' s delta defined to be u n i t y if ra=0 a n d zero if m # 0 . Obviously a new r e s i d u a l vector r ( n ) satisfies the r e q u i r e m e n t of m u t u a l i n d e p e n d e n c e for the c o m p o n e n t s since
( r ( n ) r ( n ) T) = p - t ( r ( n ) r ( n ) T ) p - T = p-1Ep -T
6) Another
f e a t u r e c a n be seen f r o m
(r(n)Tr(n))
=
the r e s i d u a l
(r(n)TE-'r(n)),
sum
of s q u a r e s 7)
w h e r e E -1 is served as a m a t r i x of w e i g h t s . Hence the above new m o d e l i n g is based on the w e i g h t e d l e a s t - s q u a r e s f i t t i n g of the generalized M A R model in E q . ( 3 ) , while an o r d i n a r y M A R m o d e l i n g is the case of a unit m a t r i x of weights a n d A ( 0 ) - E . The p r e s e n t p a p e r is o r g a n i z e d as f o l l o w s . C h a p t e r 2 is devoted to a t h e o r e t i c a l d e r i v a t i o n of the M A R model in E q . ( 3 ) as a g e n e r a l r e p r e s e n t a t i o n of a vector t i m e series. C h a p t e r 3 describes a m e t h o d a n d p r o c e d u r e for the modeling. It is s h o w n t h a t the p a r a m e t e r s A(m)(11~=O,1,2...M) a n d E c a n be d e t e r m i n e d r e c u r s i v e l y f o r e a c h value of /~1 by using c o n v e n t i o n a l a l g o r i s m s such as the LWR o n e . C h a p t e r 4 is concerned with s p e c t r a l p r o p e r t i e s of the e s t i m a t e s in t e r m s of c l o s e d - a n d o p e n - l o o p c h a r a c t e r i s t i c s . In p a r t i c u l a r ,
Multivariate autoregressivemodeling
699
it is shown that the magnitudes of residual variances can be chosen freely by adjusting the weight m a t r i x without a n y change in spectral estimates. Fi nal l y, c h ap ter 5 gives some concluding r e m a r k s .
2.
MAR REPRESENTATION
For a multivariable system to be linear and s t a t i o n a r y , consider a vector time series { x ( n ) } as a set of observations made at equidistant time intervals n h ( n = O , 1 , 2 . . . . ). In order to derive an MAR representation of { x ( n ) } , we first define a p o w e r - s p e c t r a l - d e n s i t y ( P S D ) m a t r i x P ( f ) for { x ( n ) } :
p(f)
(8)
zx ( X ( f ) X ( f ) * ) ,
where f is a f r e q u e n c y , X ( f ) is a Fourier component vector for { x ( n ) } , and the superscript t denotes the operation of H erm i t i an transposition. Since P ( f ) is H e r m i t i a n , it can be factorized into the diagonal m a t r i x E ( f ) of real positive eigenvalues and the o r t h o g o n a l m a t r i x U ( f ) of the eigenvectors:
P(f) We m a y rewrite Eq.
= U(f)E(f)U(f)' (9)
P(f)
(9)
as follows.
= U(f)E(f)~/2E(fN)-~/2E(fbl)E(fN)-~/zE(f)x/2U(f) ' = T ( f , f N ) - I E ( f N ) [ T ( f , f N ) - x ] *,
(10)
where fN is the Nyquist frequency equal to 1/2h and
T(f,fN)
= [U(f)E(f)I/2E(fN)-X/2} -I
(11)
Since T ( f , fN)-I is causal and asym pt ot i cal l y stable, it satisfies the condition of invertibility for a moving average process ( H a y k i n and Kesler, 1979). Hence we can express T ( f , f N ) in the form of an infinite Fourier series :
T(f,fN)
= ~-~,A(m)z -'h,
(12)
m=O
where z=exp(i27zf) given by
and A(m) is the Fourier coefficient m a t r i x of the expansion
A(m)
= f_tNT ( f fN
fN)z'hdf '
(m=0
1 2 ...). '
(13)
'
Hence we have expressed P ( f ) in terms of A(ra) and E ( f N ) . This leads us to represent { x ( n ) } as an MAR process in the form of Eq. ( 3 ) where the regression m a t r i x A(la) is defined by Eq. ( 1 3 ) and the residual-noise
700
N. MORISHIMA
vector r(Tt) has noise processes to be white and m u t u a l l y independent such that
(r(n)r(n')r>
=
Eti.,~,,
(14)
with the diagonal m a t r i x of variance fN = ./-/fE(fN)df
=
E(fN)/h.
(15)
Notice t h a t , taking the Fourier t r a n s f o r m of Eq. ( 3 ) and inserting the result on X ( f ) into Eq. ( 8 ) , we then obtain P ( f ) in Eq. ( 1 0 ) . Obviously the MAR expression thus derived is characterized by the z e r o - l a g regression m a t r i x A ( 0 ) and the diagonal m a t r i x P. of residual noise. This is in co n tr as t to an o r d i n a r y MAR model since the latter can be derived from Eq. ( 3 ) as follows. By s ubt r a c t i ng x ( a ) from the both sides of Eq. (3 and multiplying the result by A ( 0 ) from the left, we then obtain Eq. (1 where A(m)
= A(0)-la(m)
Y(n)
= A(0)-~r(n).
(m=l,2 ....
),
(16
and (17
Consequently the o r d i n a r y MAR modeling does not distinguish a z e r o - l a g correlation f r o m the others. This is one of the reasons why we often encounter an undesirable significant c or r e l a t i o n of residual noises in spite of a considerably short sampling period.
3.
MAR MODELING
The Y u l e - W a l k e r equation for the present MAR modeling can be derived as follows. Using Eq. ( 1 3 ) for A(m) and also the covariance function m a t r i x R ( l ) given by
R(l)
we m a y calculate Kesler, 1 9 7 9 )
the
A(m)R(l-,~
=
f_/~ S~ P ( f ) z -th d f
convolution sum
(l=0
of A(la)
'
-el and
'
+2
'
R(l)
= (i27rh)-'E(fN)f{T(f,fN)-'}*z'-ldz.
...), as ( H a y k i n
(18) and
(19)
•=0
Since {T(f,fN)-~} * is anal yt i c on and inside the unit circle in the z plane, we can p e r f o r m the integral in E q . ( 1 9 ) and arrive at the Yule-Walker
Multivariateautoregressivemodeling
701
equation
~, A ( m ) R ( l - m ) • =o
-
{ E{A(0)T]-~ 0
(l=O), (l>1).
(20)
Once the values of R ( l ) are given, all the model parameters A(m) and E can be determined uniquely by using Eq. ( 2 0 ) . To show this explicitly, we m a y rewrite Eq. ( 2 0 ) for an order M in the form n An(O)R(I)An(O)' - -~,An(m)An(O)-1"An(O)R(t-/n)An(O)' m=l
(l=1,2
....
(21)
~/),
and n
EM-
An(O)R(O)An(O)'
+ ~,An(m)An(O)-t.AM(O)R(rn)'AM(O) '
(22)
Note that these set of equations are similar in form to the o r d i n a r y Yule-Walker equations, provided that the values of AM(O)R(I)AM(O)* are given to solve An(re)AM(O) -1 and En. That is, once an adequate initial guess for An(0) and a set of R ( l ) ( l = O , 1 , 2 . . . M ) are given, we can calculate all the values of A M ( m ) ( m = I , 2 , . . M ) and EM. Such a calculation will be repeated for various values of AM(0) until a diagonal form of EM is obtained. This means that an iterative least-squares fitting is necessary for unbiased estimation o£ the model p a r a m e t e r s . However, using the following procedures, we can determine all the model parameters at once, not iteratively. ( 1 ) Define the o r d i n a r y Yule-Walker equations which are derived from Eqs. ( 2 1 ) and ( 2 2 ) : M
R(l)
= -)--~,An(m)R(/-m)
( / = 1 , 2 . . . . M),
(23)
m=l
and M
~M = R ( 0 )
+ ~--~,AM(m)R(-ra),
(24)
lI~=1
where
AM(m) = AM( O)-tAM(m) ,
(25)
EM = AM( 0 )-IEM( An( 0 )-1 ) ,
(26)
and
( 2 ) Solve Eqs. ( 2 3 ) and (24) to determine the values A M ( m ) ( m = l , 2 . . . . ~q) and EM. ( 3 ) Factorize En into eigenvalues m a t r i x D and eigenvectors m a t r i x U:
EM = UDU'. (4) obtain
Notice
the
correspondence
of
(27) between
Eqs.
(26)
and
(27),
we thus
702
N. MORISHIMA
[:M = D and also,
f r o m Eq.
and
AM(O) = U ' ,
(28)
(25),
aM(m)
ffi A M ( 0 ) ~ M ( m ) ,
(m=I,2
....
H).
(293
There is nothing special for the p a r a m e t e r calculation, since an o r d i n a r y LWR a l g o r i t h m is utilized as it is. Only a care should be taken of the factorization of EM. T hat is, the m a t r i x U must be determined in correspondence to the structure of EM. For e x a m p l e , in the case of the correlated residuals where the cross correlation~ E M ( i , j ) and E M ( j , i ) are significant, U must have eigenvector components in the elements U ( i , i ) , U(i,j), U ( j , i ) and U ( j , j ) . Such a m a t r i x structure of U can be constructed easily and uniquely, for e x a m p l e , by the use of J a c o b i ' s method for o r th o g o n al t r a n s f o r m a t i o n of a real s y m m e t r i c m a t r i x .
4.
PROPERTIES
OF
SPECTRAL
ESTIMATES
4. I. Closed-Loop Properf~ The present modeling is equivalent to the o r d i n a r y one for spectral estimation of o v e r a l l ( c l o s e d - l o o p ) p r o p e r t y . This is ascertainable in terms of a PSD m a t r i x . The Fourier t r a n s f o r m of Eq. ( 3 ) yields
P(f)
~ (X(f)X(f)¢) M
M
- ( Yq,AM(m)
) -'r.M ( EAM( m)'
M
=
)
~=0
I~=0
M
(E+E]M(m)z-'h) -' ~M (E+Y]~M'(m)z'h) -' • =1
Obviously the last result in Eq. MAR modeling.
(30)
M=I
(30)
is just the PSD given by the o r d i n a r y
4 . 2 . Open-Loop Property Distinct difference with r e g a r d to the PSD is the power contributions from each residual noise. This can be seen f r o m Eq. ( 3 0 ) in terms of the a u t o - P S D s due to each diagonal of F.M for the present modeling, and the those f r o m each diagonal of ~M for the o r d i n a r y modeling. An i m p o r t a n t point is that the f o r m e r preserves the total power of a signal by the sum of power contributions, in contrast with the latter neglecting all the power contributions f r o m the o f f - d i a g o n a l s of ~M. The other o p e n - l o o p characteristics such as a direct response f r o m one signal to an o t he r and a noise source PSD are also different between the present and o r d i n a r y MAP~ modeling. This will be shown in the following. Taking
Multivariate autorcgressivemodeling
703
the Fourier t r a n s f o r m of Eq. ( 3 ) and splitting the regression m a t r i x diagonal one A d ( f ) and an o f f - d l a g o n a l one A o ( f ) as
-EA/~(m)z -'h -
Ad(f)
into a
(31)
+ Ao(f),
m=0
we thus obtain the i n p u t - o u t p u t relation for X ( f ) : X(f)
where G ( f )
G(f)X(f)
(32)
+ F(f),
is a f r e q u e n c y - r e s p o n s e function m a t r i x given by C(f)
and F ( f )
-
-
-Ad(f)-*Ao(f),
(3 3)
is a noise-source vector defined by the a u t o - P S D m a t r i x O(f)
-
(F(f)F(f)'}
-
Ad(f)-lEM(Ad(f)-') ' .
Consequently, G ( f ) and Q ( f ) the values of A~/(ra) and EM.
4.3. Residual-Noise
(34)
depend on the method of MAR modeling t hrough
Maonitudes
As shown in Eq. ( 6 ) , the new MAR modeling yields a unit m a t r i x of residual, while the corresponding one in Chapter 3 uses the ei genm at ri x D as residual m a t r i x E,~. This suggests t ha t the magnitudes of residual noise m a y be chosen freely without a n y change in spectral estimates. To prove this, we will factorize I:• in the f o r m , instead of Eq. ( 2 7 ) ,
(35)
EM ffi UV*/~WV'nU ' ,
where V and W are under the condition
the diagonal
matrices
with
arbitrarily-chosen
D = VIn~fV I n .
(36)
T hen , f r o m the correspondence between Eqs. new p a r a m e t e r s as F.M = W so that the other same way as Eq.
and
diagonals
(26)
and
(35),
7~M( O ) = V-~nU * ,
parameters A(m)(m=l,2 .... (29). Hence we obtain
/q)
we m a y define
(37)
can
be determined in the
H
E A ~ / ( m ) z -'h = V - ' h A d ( f ) m=0
+ V-~nAo(f),
(38)
704
N. MORISHIMA
Consequently the frequency-response function m a t r i x G ( f ) PSD m a t r i x Q ( f ) are given by G(f)
and the noise-source
= -(V-~/2Ad(f ) )-~(V-t/2Ao(f) ) -
-Ad(f)-~Ao(f),
(39)
and Q(f)
= ( V - t / 2 A d ( f ) )-tlc( ( V - ~ n A d ( f ) ) - ~ ) ' -
Ad(f)-~O(Ad(f)-~) '
(40)
Obviously ~ ( f ) and Q ( f ) are the results given by Eqs. ( 3 3 ) and ( 3 4 ) , respectively, and therefore independent of the factorizatlon of ~t~ into a diagonal m a t r i x E• and an orthogonal one A ~ ( 0 ) .
5.
CONCLUDING
REMARKS
For a linear s t a t i o n a r y vector time series, the MAR model has been derived theoretically in the generalized form with a zero-lag regression matrix and independent residual noises. It has been shown that the MAR modeling can be performed easily using the o r d i n a r y algorithm for calculating recurslvely a set of model parameters up to a specified order. The present modeling would be useful in estimating open-loop characteristics, notably when a residual-nolse correlation is significant. For example, the combined use of the present and o r d i n a r y modeling would be promising for checking the effect of residual-noise correlation on open-loop estimates. Such a spectral analysis is in progress now and the results will be reported somewhere in near future.
REFERENCES
H a y k i n , S. and S. Kesler (1979) Predictlon-Error Filtering and M a x i m u m - E n t r o p y Spectral Estimation, in N o t t l i n e a r ~ e t h o d s o f S p e c t r a l A a a l ! l s i s , edited by S. H a y k i n . S p r i n g e r - V e r l a g : Berlin, pp. 9-72. Kanemoto, S, Y. Andoh, F. Y a m a m o t o , K. Kitamoto and K . Nunome ( 1 9 8 2 ) Identification of Pressure Control System Dynamics in BWR Plant by Multivariate Autoregressive Modeling Technique. J. Nucl. Sci. Technol. 19, 5 8 - 6 8 . Kleiss, E. B. ( 1 9 8 2 ) Determination of BWR Incore Power Feedback Effects by Radial Coherence Measurements. P r o q . N u c l . E n e r g y . 9, 6 6 5 - 6 7 4 . Morishlma, N ( 1 9 8 6 ) On the Frequency Resolution of Conventional and Autoregressive Spectral Estimates of Stationary Reactor Noise. Ana. Nucl. En. 13, 125-130. Wiggins, R. A. and E. A. Robinson (1965) Recursive Solution to the Multichannel Filtering Problem. J. Geophy. Res. , 70, 1885-1891 .