Copyright © IF AC Control Science and Technology (8th Triennial World Congress) Kyoto. Japan . 1981
MODELLING
STATISTICAL PROPERTIES OF MULTIVARIATE AUTO REGRESSIVE SPECTRAL ANALYSIS H. Sakai and H. Tokumaru Department of Applied Mathematics and Physics, Faculty of Engineering, Kyoto University, Kyoto 606, Japan
Abstract. This paper investigates the statistical properties of multivariate autoregressive (AR) spectral analysis. Under the assumption that the data are taken from a pure multivariable AR process, the asymptotic error covariances between the elements of the es timated AR coefficient matrices and residual covariance matrix are derived by using the periodogram method. Next, a similar analysis is performed to obtain the error covariances beth'ecn the e lements of the estimated spectral densit y matrix. If the order of the fitted autoregression is much greater than the true order, the covariances are similar to those of the multivariate version of the classical BlackmanTukey method. Keywords. Autoregressive process; Qultivariate systems; identification; spec tral analysis; signal processing; statistics. MULTIVARIATE AR PROCESS
I NTRODUCTIO N In recent years, multivariate AR spectral analysis method has been successfully applied to man y actual prob lems (Akaike and Nakagawa, 1972). But as for the statistical properties of this method, as far as the authors are aware, there are little substantial wor ks except for a comment of Parzen (1970). The study of this method is just a multivariate version of that of Akaike (1969), but for this purpose we must fully know the statistical properties of the parameter estimates of multivariate AR processes. The well known paper of Mann and Wald (1943) discussed these problems where they derived the asymptotic distributions of the estimates of the coefficient matrices of an AR process, but did not mention any of the statistical properties of the estimate of a residual covariace matrix, another important quantity in AR processes. In this paper, using the periodogram method (Sakai and others, 1979), we first derive several properties of the parameter estimates. Then using these results, we derive the error covariances of the estimated spectral density, thus generalizing the work of Akaike (1969). We also see that if the fitted order is taken as much greater than the true order, the resulting statistical properties of the multivariate AR spectral estimate resemble with those of the multivariate BlackmanTukey method. However, this result is rather hard to AR spectral analysis method, since the order can be satisfactorily estimated by Akaike's information criterion (AIC). In this respect, we briefly mention the recent result of Sakai (1981) .
Let us consider the following d variates mth order Gaussian AR process (1)
where AI' A2"'" Am are dxd AR coefficient matrices, {Ut } is a sequence of Gaussian white noise dvectors with mean 0 and covariance matrix W which from now on '~' e call as "residual covariance". That is, E[l}tl = g, E [l}tl}Il = ]V Ot, s. (\'I is positive definite.) The coefficient matrices safisfy the following YuleWalker equation [AI' A2 ,· .. , Am l
Ra
RI
R_l
Ra
R mI R m2
R R m+l m+2
Ra
= [RI' R2 ,···, Rml
(2)
where we define Rk = E [h+k~il. \Ve denote the dmxdm grand matrix in the left hand side of (2) by B whose (i,k)th block matrix is Rk_~ Similarly, the (i,k)th block matrix of p = R I is denoted by ~ik' The residual covariance satisfies
 AR
m m
(3)
The spectral density matrix F(s) of (1) is given by F(s) T
1
where A means (A
7a7
1
21T A (s)W·A 1 T
)
T
(s)
(4 )
and A(s) is defined by
H. Sakai and H. Tokumaru
708 m
jks
I Ak e(5) k=l where 1 is the dxd uni t matrix. To ensure the stationarity of (1) , we assume that A(s) = I 
m
 I
I1
k Ak z I
k=l
N
(6)
When a set of data {~l' .. · '~N} is available, we first estimate Rk by the following estimate 1 Nk T A AT Rk = ~ t~l ~t+k~t' R_k = Rk (k ~ 0).(7)
E[J (s)J ,(s)J (t)J ,(t)] n
f
p
n
2
nn
,(s)f
p
2
pp
,(t)+(2rr) f
np
(s)f, ,(s)x n p
(8)
IDN(s+t)12 +(2rr)2 f ,(s)f, (S)IDN (St)!2. . np n p I (14)
Then we have F ( s ) e j kS d s, RAk
jts
fpq(s) ~ (F(S))pq' The. second t~rm ~f (13) lS un1form 1n s,t so that 1tS contr1but1ons to the integrals below can be neglected. From the assumption that {~t} is Gaussian, the property of the fourth moment of Gaussian random variables, and (13), we obtain
(2rrN) f
1 N . N T . IN(s) = 2[ I X e Jts] [ I X e Jts ]. N t=lt t=lt
(IT J_ rr
e
t=l
Defining the periodogram matrix by
Rk =
I
DN(S) = ~nd
for I z I ~l.
+0
where D(.) is the Dirichlet kernel and is defined by
rr jks IN(s)e ds. rr (9)
The contribution of the first term to (12) is zero by (2). Noting that IDN(s)1 2/2rrN .... o(s), the contribution of the second term is
THE ERROR COVARIANCES OF THE ESTIMATED AR COEFFICIENT MATRICES Subst~tuting Rk into Rk in (2) gives the estimate Ai of Ai. By the same method as in Sakai ~nd others (1979), the estimation error~ 6Ai= AiAi(i=l, ... ,m) are asymptotically expressed as
a. (s)f (s)a. (s) 11n np 12 P
m
m
I
I
AkR l _k ,···, R k=l m
k=l
Ak R k]· m Noting F(s) = FT(_s), B~k' maining part of (15) becomes
From (5), (9) we have m
6Ai 
But n,p
[M l ,··· ,Mm] R
[RI 
j xa. (s)f" (s)(P k , ),. e (kk')sds. 12 P n p  m2 p )2 (15 )
f~rrA(s)IN(s) I ~ki k=l
e
jks
(10)
ds.
I
k,k',n',p'
We denote the (i,k)th element of A(s) by (A(s))ik = aik(s). We rewrite (8) as
I
k' ,pt
Ik
(~km
1
~k'k' the re
)n'J' (Pk'm )p'J' CBk'k)p'n' 1 2 2
~k' k ~ km ) p , J' (~k 'm ) p , j 1
1
2
2
(17)
Then, the error covariance between the estimated coefficient matrices is given by
where !ik is the (i,k)th block matrix of dmxdm unit matrix! and we used the relation ~.p = !. By the same argument, we have the contribution of the third term as ) ,. a. (s) J (s)J ,(s) (P kml n J 11n n n l )p" e a i p ( t ) J (t)J ,(t) (P k J2 2 p P  'm 2
~j ks
I
I
k' p,p'
jk't
2 rr N
}dsdt. (12)
It is well known (Brillinger, 1975) that E [J (s)J (t)] = 2 rrf (s) DN (s +t) + 0 (1) P q pq
X(A(S)F(S)P
) . . ej(k+k')s ds.
 kml 12 J l
(1 ~ k, k' :;, m, 1 ,:;, n, n' ,p, p' ,:;, d)
(13)
(18)
Rewriting the integral (18) into a complex integral by the transformation z exp(js) and using Cauchy's theorem, we can easily see that this integral is zero, since (4),(S),and (6), 2rrA(s)F(s) = W.AT(s) = W·Adj(IAIZ ...
Statistical Properties of Multivariate Autoregressive Spectral Analysis Amzm)/IIAIZ ... Amzml is regular inside the unit circle and k+k'l > O. Thus we have
THE ERROR COVARIANCE OF THE ESTIMATED RESIDUAL MATRIX In this section, we derive a new result about the statistical properties of the estimated residual matrix. Since the estimate is given by (20)
709
(26) This means that the estimation errors of the AR coefficient matrices and the estimation error of the residual matrix are asymptotically uncorrelated. This is a natural generalization of the corresponding result of Mann and Wald (1943). Substituting (10) into flA. in (21) and using (25), fiN is expressed as 1 T!
f _T!A(s)IN(s)A
flW·
T
(s)ds  W.
(27)
we have
Thus
the estimation error is asymptotically expressed as m
L (flA 1. R1.
i=l
+ A.flR . ) 11
T
·A (t)) . . ]ds·dt  E[(W) . . JeW) . . 12J 2 1 1J l 1 2J 2
m
f:T!A(S)IN(s)ds 
L
i=l
A.R .  W. 1 1 (21)
(W). l
Hence, we have
. E [ (W) . . ] + (N) .
lJl
1 2J 2
1
. (\'I) . . .
1J l
1 2J 2
(28)
Writing the first term of (28) componentwise, using (14), and performing the same calculation as in (12), we have
(M
) . . (M.R .) . . ]  E[(M ) . . JeW) . . • 1 1 12J 2 ml l l J l ml 11J l 1 J 2 2 This is also a generalization of the classical result of Mann and Wald. For more general (22) treatment, we refer to Ljung and Caines(1979) . Substituting (10) into the first term of (22) and using the same technique as in the calculation of (21), the first term is asymptotiTHE ERROR COVARIANCES OF MULTIcally equal to VARITE AR SPECTRAL ESTIMATE Here we examine the error covariances of the consistent estimate of (4) 1
~
~l
F(s) = f i A
From (19), the second term becomes
~
~T
(s)N'A
(s)
(30)
~
where A(s) is defined by replacing Ak in (5) by Ak' Invoking a general theory of Mann and \'laId (1944) mentioned by Akaike (1969), and noting the consistenc~es of Ak, Wk, the estimation error flF(s) = F(s)F(s) is asymptotically expressed as
L
i,p
flF(s) • A +(W). . (AT) j 1 j 2
1112
ml
1 1 T (s)M(s)F(s)+ f i A (s)6\~·A (s)
(24) F(s)flAT(s)A T (s)
where we used the relation RT. 1 ~ki' and
(31)
where m
M(s) (25)
Hence, we obtain
1
=
L
AA Ll
k e
j ks
.
(32)
k=l Since from (26), flW and lIA(s) are uncorrelated, we have
H. Sakai and H. Tokumaru
710
then by the same calculation as in (35), (33) eventually becomes
E[( L'.F (s)) . . (L'.F(t)) . . ] l J 1 J 2 2 l l 1 1 T 1 E[2 (A (s)L'.l'i·A (5)) . . (A (t)L'.W· 4rr l lJ l AT(_t)) . . ]+E[(Al(S) L'. A(S)F(S)+F(S)L'.AT(s)' 1
2J 2
(A1(s)WAT(t)) . . (FT(s)P(s,t)F (t)) . . 1112 JlJ2
(Al(_s)\v.AT(_t)) . . (F(s)P(s,t)FT(t)) . .
AT(_s)) . . (Al(t)L'.A(t)F(t) + F(t)L'.AT(t)AT(t)) . . ]. 1 J 2 2
(33)
Using (19) and (29), we calculate (33). Putting Al(s) = C(s) = (cik(s)), the first term of (33) is
C. l
ln
1
T
T
T
(s)W ' A (t)) . . (F (s)P(s,t)F (t)) . . 1
2J 2
Jl
1z
+(A 1 (t)N·A T (5)) . . (F T (t)P(t,s)F T (5)) . . 1
J2~
2J l
1 1 T lT + 2(A (s)\v·A (t)) . . (A (s)W·A (t)) . . J 4rr l l 2 3t1z (37)
c.1 n ,(t)(L'.W), k' c J. k,(t)] n, 2 2
L
+(A
1 1 T 1 T + 2(A (s)W·A (t)) . . (A (s)l'i'A (t)). 4rr 1112 J 1 j2
_1_ E[ L C. (5) (L'.\V) k c. kCsi L n 4rr2 n,k l l n Jl n' ,k'
n,k,n' ,k'
ll1z
J lJ 2
l Jl
l
+
(s)c. kCsic . ,(tic . k,(t) 12n Jl J2
It is easy to s ee that (37) r educes to (37) of Sakai and others (1979). Using the method of Makhoul (1978) for univariate case, we easily obtain a Cholesky decomposition of B 1 as
[(Al(S)W'AT(t)) . . (Al(_s)W.AT(_t)). . 1112 JlJ2 + (A
1
(s)l'i·A
T
(t)) . . (A l J l 2
1
(s)\V'A
T
Also, one of the four terms arising from the second term of (33) becomes as follows;
C.
n,k
L
c.
n', k' ~
l
_1 N
1
2n
ln
\ i
(5) L (M . ) k e
,et)
\L
1n
jis f · (5) kJ l
L (M.,) , k' i' n, 1
e
ji's
1
r\
(t)) . . ]. IAI 1 J l 12 m(34)
= E [L\
B
l~~ mI
o
0
o
Ai m2
'T 'T Am_l ml Am2 m2' . . r AI mI A 2 mI I AI m2
..,1 "0
Am_l mll A I.2 .2
x
o
j
(38)
f ,·
k J (5)] 2
where Aik (i=l, ... ,k; k=l, ... ,ml) and l'Ik (k =0, . . . ,ml) are theoretical coefficient and residual matrices arising from kth order AR model fitting . From (38), P(s,t) becomes
(s)c. ,(tlCW) ,x 1 2n nn n, k ,n 'k , ,1ln C.
mI P(s,t)
L
ej(mk)(s+t)A~(S)W~lAk(t) (39)
k=O where k
Ak(S) = I (FT (5)
L
(35)
i, i'
Defining a matrix P(s,t) by P(s,t) =
\
P
L  ii' i, i'
ej(is+i't)
,
(36)
L
i=l
Aik e
j is
.
(40)
Now we consider the case where the fitting order M is much greater than the true order rn, i.e., M »m. By putting Ai = 0 (i=m+l, ... , M) we can regard {~t} as an Mth order AR process, so that (37) becomes valid after replacing m by M. But, ohviously, Ak(s)=A(sl,
Statistical Properties of Multivariate Autoregressive Spectral Analysis
711
II'k=W (k=m, .. . ,M) follow so that (39) becomes
o Theory • Simul at ion
P(s,t)

10
Since we assume M» m, first term of (41) is much smaller than the remaining one. Neglecting this term and substituting (41) into (37), we have
o
a
(_1_ Al(S)W'AT(t)) . . (_l_Al(_s)II" 2rr 11122rr T A (t)) . . (1 + 0l~_r.,(s+t) + 0M_m(st)) +

J 1J 2
o
(  12 Al(S)W.AT(_t)) . . (  12 AI (t)II" rr I J rr l 2 T A (s)). . (1 + OM (st) + O~I (ts)). (42) 1 J m m 2 l Noting that 0N(s) = N for s=0, ±2rr, ... , and 0(1), otherwise, we see that (42) is of order M only when s = it. From (nF(s))ik=(nF(s))ki, it is sufficient to consider the case s=t . In this case, (37) is expressed as
•o
. o
0.1
o
o Fig. 1
s
+
TI
Theoretical and experimental values of N.E[(nF(S))rl]'
2M(f . . (s)f . . (s)+f . . (s)f . . (s)) 12J l 1112 JlJ2 I J l 2 (s=O,rr) N·E[(nF(s)) . . (nF(t) . . ] .. 0M(l) (st±t). I J 1 J 2 2 l l ( 43)
TI/2
o Theory 100
Si mula t ion
•
~
where 0M(l) means that this quantity is of order 1 about M. o
The result takes the same form with that of the multivariate BlackmanTukey procedure (Hannan, 1970) and is a generalization of Sakai (1979). However, it is obvious that actual statistical variabilities are much smaller than (43), since we can obtain a reasonable fitting order near m by using FPE or AlC. In fact, it is shown in Sakai (1981) that the probability of selecting the correct order tends to one as the number of variates increases. For example, for two variates case, AlC selects the true order with probability 0.88 while in univariate case it is only about 0.71.
• o
o
• 10
... 0
0
o

.0 0
•
•
•
0
•
0
• 0
0
0
a
• ••• 0
0
SHruLATION
0
0
• To see the validities of the above results, numerical simulations were performed. 450 sets of data each of N = 1000 length were generated by the following two variates third order AR process
0
0 0
I
0
+ 5
TI/2
TI
Fig. 2 Theoretical and experimental values of N.E[(nF(s))~2]'
H. Sakai and H. Tokumaru
712
r 0.3
=
!'t
L +
0
[0.1 0.5
_0°5] tl X
: ] \3
[0.2
+
2.0 +
:]
~t2
U t
wi th W = diag (1. 0,1. 0). Figs 1 and 2 show the theoretical values (37) and simulation values. The agreements are good and justify the present asymptotic analysis. CONCLUSION This paper has presented the statistical properties of multivariate AR spectral analysis under the assumption that the process is finite order autoregressive. It has been shown that if the fitting order is much greater than the true order, the asymptotic properties resemble with those of multivariate BlackmanTukey method. It is a future problem to generalize the work of Berk (1974) which treats a certain class of infinite order scalar AR processes to multivariate case.
REFERENCES Akaike, H. (1969). Power spectrum estimation through autoregressive model fitting. Ann. Inst. Statist. Math., 21, 407419. Akaike, H., and T. Nakagawa (1972). Statisti
cal Analysis and Control of Dynamic Systems, SaiensuSha, Tokyo. (in Japanese) Berk, K.N. (1974). Consistent autoregressive spectral estimate. Ann. Statist., 2, 489502. Brillinger, D.R. (1975). Time Series: Data Analysis and Theory, Holt, Rinehart, and Winston, Inc., New York. p. 73. Hannan, E.J. (1970). Multiple Time Serie8, Wiely, New York. p. 280. Ljung, L., and P.E. Caines (1979). Asymptotic normality of prediction error estimates for approximate system models. Stocha8tics, 3, 2646. Makhoul, J. (1978). A class of allzero lattice digital filters : properties and applications. IEEE Trans. Acoust., Speech,
Signal Processing, ASSP26, 304314. Mann, E.B., and A. Wald (1943). On the statistical treatment of linear stochastic difference eqJations. Econometrica, 11, 173220. Mann, E.B., and A. \~a1d (1944). On stochastic limit and order relationship. Ann. Math. Statist., 14, 217226. Parzen, E. (1970). Multiple time series mode1ing. in Multivariate Analysis II, P.R. Krishnaiah Ed., Academic Press, New York. pp. 389410. Sakai, H. (1979) Statistical properties of AR spectral analysis. IEEE Trans. Acoust.,
Speech, Signal Processing, ASSP27, 402409.
Sakai, H. (1981). Asymptotic distribution of the order selected by AIC in mu1tivariate autoregressive model fitting. Int. J. ControZ, 33, 175180. Sakai, H., T. Soeda, and H. Tokumaru (1979) On the relation between fitting autoregression and periodogram with applications. Ann. Statist., 7, 96107.