A new algorithm for recursive estimation of ARMA parameters in reactor noise analysis

Ann. nucl. Energy, Vol. 19, No. 5, pp. 287-301, 1992 Printed in Great Britain

0306-4549/92 $5.00+0.00 Pergamon Press plc

A NEW ALGORITHM FOR RECURSIVE ESTIMATION OF ARMA PARAMETERS IN REACTOR NOISE ANALYSIS TRAN DINH TRI CRIP Institute for Atomic Energy Research Applied Reactor Physics Department H-1525 Budapest 114, P.O.B. 49, Hungary

( Received 20 January 1992 ) Abstract - - In this paper a new recursive algorithm for estimating the parameters of the ARMA model from measured data is presented. The Yule-Walker equations for the case of the ARMA model are derived from the ARMA equation with innovation. The recursive algorithm is based on choosing appropriate form of the operator functions and suitable representation of the (n + 1)-th order operator functions according to ones with lower order. Two cases, when the order of the AR part is equal to that of the MA part and the general case, were considered.

1. INTRODUCTION The ARMA (Autoregressive Moving Average) model can be derived from the physical model (the linear Langevin equation). Moreover for the system in a steady state and in a normal case it can be carried out either by the physical way or by the mathematical way (Kishida, 1982). In the physical way, the physical state equation, which is obtained from the linear Langevin equation by coarse graining in time, first is transformed into the non-Markovian equation by elimination of irrelevant variables. Then by projecting the non-Markovian equation on the observable variable space we obtain the ARMA model. In the mathematical way, the order of manipulations is inverted. First by projecting the physical state equation the Kalman filter will be obtained, after that by contraction of information the ARMA model will be derived from the Kalman filter. The ARMA model, which is derived from the physical model in the physical way on the elimination path, is given in the following form: n

p

Aoy(k ) + E A i y ( k - i ) = B__of(k) + E B i f ( k - i) i=1

y(k) is f(k) is a

where

(1)

i=1

a q-dimensional observable vector, q-dimensional random noise vector,

n and p are the order of the AR and MA parts, respectively, A i and B i are parameters of the AR and MA parts, respectively. In this ARMA model the physical random noise force f(k) and its statistical quantities cannot be identified fl'om measured data. This difficulty can be overcame by using an 287

288

TRAN DINH TR1

equivalent quantity called the innovation in the time-series analysis. This innovation is statistically equivalent to the physical random noise in the sense that both these quantities generate the same output y(k), and that from them the same correlation functions and power spectral densities will be obtained (Kishida, 1984). So, instead of Eq. (1) we use the following ARMA model with the innovation, which has been obtained by projecting the non-Maxkovian representation (1) on the observable variable space as mentioned above: p

y(k) =

Aiy(k

-

i) + 7(k) + E BiT(k

i=1

-

i)

(2)

i=1

where 7(k) is the innovation: V(k)=y(k)

y(klk-1)

-

(3)

y(k I k - 1) is the least-squares estimator, i.e. the most probable estimator under the condition that 1) = [yT(k 1),yT(k-- 2),...] (4) y T ( k

-

-

is given, T denotes transposition. For the system in a normal case, the state variables possess the Gaussianity (Tomita, 1974; Kishida S. et al., 1976). Then the value of estimator y(k I k - 1) can be evaluated by the Gaussian conditional value:

u(k I k - i) = < y(k) I Y(k - 1) > = < y ( k ) Y r ( k - I) > < Y ( k - 1 ) Y r ( k - I) > - 1 Y ( k - I)

(5)

where < ... > is the ensemble average. The parameters of the AR model, that is a special case of the ARMA model when P

we truncate the MA term y] B i f ( k - i) in the ARMA Eq. (2), can be determined fi'om i=l

measured data by the Levinson-Durbin algorithm (Fukunishi, 1977; Upadhyaya et al., 1980; Whittle, 1963). However, in reactor noise analysis for the ARMA model there is no similm' practically applicable algorithm. In this paper we propose a new recursive algorithm for practically estimating the parameters of the ARMA model from measured data. In section 2 the system of the Yule-WMker equations for the parameters of the ARMA model will be derived. Section 3 deals with the new recursive algorithm. In section 4 we introduce an important theorem, based on which we can show the meaning of the innovation in Eq. (3) and the relation between the parameters of the AR and ARMA models. In the final section 5 some concluding remarks will be madel

Recursive estimation of ARMA parameters

289

2. S Y S T E M OF T H E Y U L E - W A L K E R EQUATIONS F O R T H E P A R A M E T E R S OF

T H E ARMA M O D E L In this section we give the system of the Yule-Walker equations, which allows us to determine the parameters of the ARMA model. Introducing the time shift operator z -1 and z as: z-ly(k)=y(k-1)

;

zy(k)=y(k+l)

,

(6)

we can rewrite Eq. (2) in the following form:

A(z)y(k) = B(z)7(k ) ,

(7)

where the operator functions A ( z ) and B ( z ) are defined as:

A(z)= I-

A(i)z -i ; i=1

B(z) = I + E B(i)z-i ' i-----1

(8)

and I is the q-dimensionM unit matrix. Using a matrix M defined as: M=
>
>-'

,

(9)

the estimator y(k t k - 1) in Eq. (5) can be rewritten as:

y ( k l k - 1) = M Y ( k -

1) .

(10)

Substituting Eq. (10) into the definition (3) of 7(k) we have:

7(k) = y ( k ) - M Y ( k - 1) .

(11)

To find equations for parameters A(i) and B(i) we begin with taking the ensembh ~ average < . . . . y T ( k -- j ) > of b o t h sides of Eq. (7). T h e n we have:

A ( z ) < y(k)yT(k - j )

> = B ( z ) < '7(k)yT(k - - j ) >

(12)

Substituting the innovation (11) into this equation we obtain that

A ( z ) C ( j ) = B ( z ) [ C ( j ) - M < Y ( k - 1)yT(k - j) >] ,

(13)

where C ( j ) is the correlation matrix:

C(j) =< y(k)yT(k - j) .

(14)

290

TRAN DINH TRI

Let us consider the term M < Y ( k - 1)yT(k -- j) >. We can see that

u(k- 1) } < Y(k-

1)yT(k - j) > = < {

y(k .- 2)

{

c(j - 1) / c (j .- 2)

YT( k _ J) > =

y(k - ' N a n )

(15)

/

C(j -" WAR)

where NAR is a positive integer, which was used to determine the innovation 7(k). The meaning of NAR and the choice of it will be discussed in section 4. It is easy to see that M = < y ( k ) y T ( k - 1) > < Y ( k - 1)yT(k - 1) >-1

= [C(1),C(2),...,C(NAR)]

C(0) C(-1) C(-2)

C(1) C(0) C(-1)

•

...

C(NAR

C(1) C(2)

... ...

C(NAR -- 2) C(NAR--3)

.

C(--NAR + 1) = [Ma,M2,...,MNas]

C(2)

.

. .

......

•

C(-1)

(~'(0)

,

)

-- 1 ) \

-1

(16)

where Mi are (q x q) matrix determined by the correlation matrices C(j). From (15) and (16) we have: CI(j-1) M < Y ( k - 1)yT(k --j) > = [M1, M 2 , . . . , MNAn]

}

6 ( j .-- 2)

=

C(j - WAR) NAR

NAR

M i C ( j - i1 = E

= E i=l

Miz-iC(J) "

(17)

i=1

Now substituting the expression (17) into Eq. (13) and changing j to k we obtain the following equations: NA R

A(z) - B(z)(Z - E

Miz-i)

C(k) :

0

,

(1S)

i=l

where k = 1,2,3,... The set of the equations (18) is the system of the Yule-Walker equations for the A R M A model.

Recursive estimation of ARMA parameters

291

3. R E C U R S I V E A L G O R I T H M F O R T H E P A R A M E T E R S O F T H E ARMA M O D E L In this section we construct a new recursive algorithm for solving the Yule-Walker equations (18). We will consider this problem for two cases: p = n and p = n - r, where 1 < r < n and n and p are the numbers of the AR and MA terms, respectively, according to EQ. (2). 3.1. T h e c a s e p = n T h e solution of the Yule-Walker equations (18) will be found by a recursive procedure. So, first let n and p be equal to 1. T h e n the Yule-Walker equations become:

(I-

Al(1)z -1) - ( I + B l ( 1 ) z - a ) ( I -

Na_a ] Z M i z - i ) C(k) = 0 , i=1

(19)

where the subscripts of A~(1) and Bi(1) indicate that n = I. To obtain A~(1) and BI(1), we put k = 1 and k = 2 into Eq. (19) and solve for AI(1) and BI(1), then we have the results: NAN

NAR

WAR

Bx(1) = [~-~ M i C ( 2 - i ) C - ' ( 1 ) - Z M i C ( 1 - i)C-1(0)][~--~ M i C ( - i ) C - I ( O ) i=1 i=1 i=1 NAR -

Z

MiC(1 - i)C-!(1)1-1

(20)

i=1

NAR

NAR

AI(1) = Z MIC(2 - i ) C - a ( 1 ) - BI(1)[I - ~ M i C ( 1 - i ) C - l ( 1 ) l . i=1 i=1

(2i)

For n = 2, the second order operator functions A2(z) and B2(z) will be seeked in the following form: A2(z) = Aa(z) + Da-Al(z)z -2 (22.a) B2(z) = S l ( z ) + D~-B~(z)z -2 ,

(22.b)

where D~ is a constant matrix and D1 # 0, Al(Z) and B~(z) are the so called adjoint operator functions: Al(z)= I-

Al(1)z ;

B,(z) = I + B,(1)z .

As seen later, these functions satisfy the following equations: NAR

{ [ I - A l ( 1 ) z ] - [I + B l ( 1 ) z ] [ I -

~ Miz-i]}C(k) = 0 , /=1

k--o & -]..

(23)

292

TRAN DINH TRI

By solving these equations for A~(1) and BI(1) we obtain: NAR

NAR

NAr~

B~(1) = [~--~ M~C(-1 - i)C-~(O) - ~ i= 1

MiC(-i)C-~(1)I[~-~ MiC(1 - i)C-1(1) -

i=1

i=l

NAR

-

~ M,C(-i)C-I(O)1-1

(24)

i=1

and NA R

NA R

M i C ( - i ) C - I ( 1 ) - B~(1)[I - E

A~(1) = E i=1

MiC(1 - i)C-1(1)1 •

(25)

i =1

Now assume that the upto n-th order solutions of the Yule-Walker equations (18) have been determined. Then the n-th order functions An(z) and Bn(z), which are defined by: n

n

An(z) = I - E An(i)z -i ," Bn(z) = I + E Bn(i)z-i , i=1

(26)

i=1

are known and they satisfy the following n-th order equations: NAR

]

An(z) - B , , ( z ) ( I - E,=, Miz-i)] C(k) = O ,

k=1,2,...,n+l

,

(27)

where An(i) and Bn(i) are the n-th order coefficient matrices (i.e. parameters of model). The (n + 1)-th order operator functions An+a(z) and Bn+l(Z) will be determined according to the n-th order functions in the following form:

A,+l(z) = An(z) + DnA,(z)z -
= B.(z)

+ D.-~.(z)z

-("+1)

(28.a) ,

(2S.~)

where Dn is an unknown constant matrix and An(z) and Bn(z) are the n-th order adjoint operator functions. To find equations, which determine An(z) and Bn(z), we substitute the expressions (281) into the (n + 1)-th order equation, then we have: NA R

JAn(z) - Bn(z)(I - E

Miz-i)]C(k)+

i=1 NAR

+Dn[A,(z) - "Bn(z)(I - E i=1

Miz-i)]z-(n+l)C(k) = 0 .

(29)

Recursive estimation of A R M A parameters

293

We can see that in Eq. (29) the first term for k = 1 , 2 , . . . , n + 1 is equal to zero because this term is the left side of the n-th order equation (27). Hence NAR

D n [ A , ( z ) - -Bn(z)(I - E

Miz-i)]z-(n+a)C(k) = 0 ,

k= 1,2,...,n+1

(30)

.

i=1

Since D . # 0, with respect to the definition (6) of the oprator z -i, we have: NA It

JAn(z) - -B,(z)(I - E

Miz-i)]C(k - n - 1) = 0 ,

k = 1,2,...,n + 1 .

(31)

i=1

Accordingly

] J A n ( z ) - -B,~(z)(I - L M i z - i ) |

C(k)= 0 ,

k = 0,-1,... ,-n

.

(32)

]

i=1

These equations are called the n-th order adjoint equations. The adjoint functions are defined as: n

-£(z)

= I -

i ;

= I +

i=l

-gn(i)z i ,

(33)

i=1

where A,(i) and B , ( i ) are the adjoint coefficient matrices. Eq. (29) for k = 1 , 2 , . . . , n + 1 is satisfied for arbitrary Dn. However, it is not satisfied for k = n + 2, thus Dn will be determined by putting k = n + 2 into Eq. (29) to obtain: NAR

D , = - {JAn(z) - B , ( z ) ( I

- E

M i z - ' ) l C ( n + 2)}x

NAR

x {[A,(z) - Bn(z)(I - E

Miz-illC(1)}-I

(34)

i=l

T h e (n + 1)-th order adjoint operator functions will be determined by expressions similar to the expressions (28): A . + l ( Z ) = A . ( z ) + E . A n ( z ) z n+l (35.a)

-B.+I(z) = -B.(z) + E . B n ( z ) z n+l ,

(35.b)

where En is a constant matrix to be determined later. Note that An+l(Z) and B . + l ( z ) are the adjoint functions, they consist of different powers of the forward shift time operator z, so in the expressions (35) we have used the factor z n+l instead of z -(n+l). It is worthy to mention that the form of the (n + 1)-th order operator functions given in expressions (28) and (35) guarantees that if we put them to the (n + 1)-th order equations

294

TRAN DINH TRI

(27) or (32), then the left side of each of these equations breaks up into two n-th order terms. This was shown in EQ. (29) for the form of A,+a(z) and B , + l ( Z ) . Now it will be shown for A , + l ( Z ) and Bn+a(z). In fact, putting the expressions (35) into the (n + 1)-th order adjoint equation we have:

NAR -

-

M,z-')lC(k)+

i=1

NAR +E,[A,(z) - Bn(z)(I - E

Miz-i)]zn+lC(k) = 0 .

(36)

i=1

In Eq. (36) the first term for k = 0, - 1 , . . . , - n according to Eq. (32) is equal to zero. Since E , # 0 and the factor z"+aC(k) = C(k + n + 1) for k = 0 , - 1 , . . . , - n can be written as C(k) for k = 1, 2 , . . . , n + 1, we can see that the second term without E , in Eq. (36) is the left side of the n-th order equation (27) and so it vanishes too. The equation (36) for k = 0, - 1 , . . . , - n is satisfied for arbitrary E , . However, it is not so for k = - n - 1 , then the constant matrix E , can be obtained by putting k = - n - 1 into Eq. (36):

NAR E, =-

{[A,(z)--B,(z)(I-

E

M i z - i ) ] C ( - n - 1)} x

i=1

NAR

× {[an(z) - B,(z)(I - E

Miz-i)lC(O)}-a

(37)

i=l

To get recursive formula for the coefficient matrices A,+l(i), we substitute the expression (27) and (33) into the expression (28.a), then we have: n+l

n

n

I-- E An+l(i)z-i = I - Z An(i)z-i + D , [ I - E An(i)zilz-("+U i=1

i=1

i=1

By comparing coefficients of each z -i in both sides of the above equality we obtain the recursive formula for An+l (i): D

An+l(i) An+l(n + 1)

=A,(i)+D,A,(n+l-i)

,

= -D.

i=1,2,.,

n " '

(38)

In the same way we obtain the following recursive formulae for B . + I ( i ) , An+l(i) and B.+I(i):

i

B,+l(i) B.+l(n + 1)

{A,,+I(i) A , + l ( n + 1)

= B,(i) + D , B , ( n + l - i) ,

i = l,2,...,n

(39)

= D.

= A,(i) + E , A , ( n + l - i) , = -E,

i= 1,2,...,n

(40)

Recursive estimation of ARMA parameters

{

=B,~(i)+E,,Bn(n+l-i),

B,.,+,(i) ~ , , + l ( n + 1)

295

i=l,2,...,n

(41)

= E.

3.2. T h e g e n e r a l c a s e p = n - r We will consider only for 0 < r < n, because r = n is the case of the A R model, and r = 0 is the considered case p = n. T h e n - t h o p e r a t o r functions axe defined as: ~

n--T

A.(z) = I-

A.(i)z -i ;

Bn(z) = I + Z

i=1

B'~(i)z-i

(42)

i=1

T h e solution of the Yule-Walker equations (18) will be found in the recursive way, then we have to determine sequentially the first, second, . . . order solutions. In addition, for n = 1 , 2 , . . . , r we take B , ( z ) = I , and we denote such values of n by n'. As shown later in section 4 NAR

[I- Z

MJ z - j l C ( k ) = O

,

k= I,2,...,NAR

•

j=l In addition, we suppose t h a t r < NAR, then for n' the Yule-Walker equations (18) can be written as: A',(z)C(k)=O , k=l,2,...,n' (43) These are the Yule-Walker equations of the A R model. So their solution can be obtained by the Levinson-Durbin algorithm (Fukunishi, 1977; U p a d h y a y a et al., 1980). T h e prime symbol at A' indicates t h a t this quantity belongs to the A R model. We continue to consider for n = r + 1, r + 2 , . . . First, for n = r + 1, the (r + 1)-th solution will be found in the following form:

A,.+,(z) = A',.(z) + Dr-A,.(z)z -(r+l)

(44.a)

Br+l(Z) = I + D,.z -1 ,

(44.b)

where D r is a constant matrix. By substituting the expressions (44) into the Yule-Walker equation (18) we obtain the following equation for A~(z): NAR [A~(z) - z ~ ( I -

Z

Miz-i)]C(k) = 0 ,

k = -1,-2,...,-r

.

(45)

i=1 To solve Eqs. (45) in the recursive way, it is required to solve the following equations: NAR

[~,,(z)- zr(I-

Z i=l

Miz-i)lC(k) = 0 ,

k = -1,-2,...,-n'

,

(46)

296

TRAN DINH TRI

wheren 1=1,2,...,r. The matrix ~1 (1) of the first order solution ~1 (z) is given by:

NAR a,(1) = [ C ( - 1 ) -

C(r - 1) + E

Uiz-iC(r - 1 ) ] C - ' ( 0 ) .

(47)

i=1

The (n + 1)-th order solution will be found according to the n-th order solution in the following form: -d,,+l(z) = "d,,(z) + Hn,an,(z)z '¢+' , (48) where H , , is a constant matrix. Substituting (48) into Eq. (46) for n' + 1 we have:

NAR [K,,(z) -

z~(I - E Miz-i)]C(k) + gn'a"'(z)zn'+'C(k) = 0

(49)

i=1

k=-l,-2,.

• .

,-n',

--~r'$

t

From Eq. (49) we can see that a,,(z) = A~,(z), where Hn, is obtained by putting k = - n ' - 1 into Eq. (49):

-1

.

A~,(z) belongs to the AR model.

NAR Hn, = -{[~n,(z) - z r ( I - E M i z - i ) ] C ( - n ' - 1)}{d~'(z)C(O)}-I

(50)

i=1

The recursive formula for coefficient matrix ~,,,+1(i) is the following: {a,v+l(i) ~ n , + l ( n ' + 1)

=-dn,(i)+H,vA~,(n'+l-i)

,

i=1,2,.,

= -Hn,

n' " '

So we have obtained the function A~(z) = ~ ( z )

(51)

.

The n-th order solution satisfies the following n-th order Yule-Waker equation:

NAR

]

A n ( z ) - B n ( z ) ( I - E Miz-i) C(k)=O ,

k = 1,2,...,n .

(52)

i=l

Now we assume that the solutions A~+l(Z) & B~+l(z), A~+2(z) & B~+2(z),..., Bn(z) are known. Next we have to find the (n + 1)-th solution of the equation

NAR

]

A n + l ( z ) - B n + , ( z ) ( I - E Uiz-i)[ C(k) = 0 , i=1

A,,(z) &

k = 1,2,...,n+

1 .

(53)

.I

The (n + 1)-th operator functions are determined in the following form:

An+l( z) = An(z) + Dn--An(z)z-('~+1)

(54.a)

Recursive estimation of A R M A parameters

297

Bn+l(z) = Bn(z) + D,~B,~(z)z -(n+l-')

(54.b)

where D , is a constant matrix.

An(z) and Bn(z) are defined by: -~,(z) = I - ~ - ~ , ( i ) z

~;

~.(z)=I+~_~n(i)z

i=1

~

•

(55)

i=1

Substituting the expressions (54) into Eq. (53) we have: NA~

[An(z)-Bn(z)(I - E

Miz-i)lC(k) + D n [ A n ( z ) -

i=1 NAR

--9.(z)z~(I-

~

M,z-i)lz-¢"+l~C(k) = O ,

k = 1,2,...,n

+ 1 .

(56)

i=1

In this equation, for k = 1, 2 , . . . , n the first t e r m of the left side is equal to zero, then the second t e r m is zero too. Since Dn # 0 and z - ( n + l ) C ( k ) : C ( k - n - 1) can be written as C(k) for k = - 1 , - 2 , . . . , - n , we have: NAR

[An(z) - B , , ( z ) z r ( I - Z

Miz-i)]C(k) = 0 ,

k = -1,-2,...,-n

.

(57)

i=1

T h a t is the n - t h order adjoint equation for the n-th order adjoint functions An(z) and Bn(z). T h e constant m a t r i x Dn is determined by putting k = n + 1 into Eq. (56): NA R

Miz-J)lC(n + 1)}x

Dn = - {[An(z) - Bn(z)(I - E j=l

NAR

x {JAn(z)- Bn(z)zr(I- Z

Miz-J)]C(O)}-I

(58)

j=l

T h e (n + 1)-th order adjoint functions are determined by

-An+,(z) = -An(z) + EnA,(z)z "+1

(59.a)

-Bn+,(z) = -gn(Z) + E , B , ( z ) z n+~-~ ,

(59.b)

where En is a constant matrix. Similar to the case p = n, if we substitute the expressions (59) into the (n + 1)-th order adjoint equation for A n ( z ) and B n ( z ) , then we have: NAR

JAn(z) - -Bn(z)z"(I - E

Miz-i)]C(k) + En[An(z)-

i=1 NAIl

- Bn(z)(Z - ~ i=1

Miz-')lzn+'c(k)

= 0,

k = -1,-2,...,-,~

-

1 .

(60)

298

TRAN DINH TR1

We can see that in this equation for k = - 1 , - 2 , . . . , - n both the first and second terms are equal to zero automatically. The constant matrix En is determined by putting k = - n - 1 into Eq. (60): NA R

E, = - {[A,(z)--B,(z)zr(I

- E

Miz-i)]C(-n

- 1)}x

i~l

NAR

x {[An(z)- B n ( z ) ( I -

E

(61)

M'z-')]C(O)}-I

i=l

To get the recursive formulae of A,+x(i), Bn+~(i), An+a(/) and Bn+~(i), we substitute the expessions (42) and (55) into the expressions (54) and (59), then by comparing coefficients of each z -i or z i in both sides we obtain the following recursive formulae:

{ {

An+l(i)

= A n ( i ) + D , - A n ( n + 1 - i)

An+l(n + 1)

=-Dn

A.+l(i) -~.+~(n + 1)

B,,+a(i) B n + l ( n + 1 - r)

,

i =

1,2,...,n

(62)

= A . ( i ) + E , , A . ( n + 1 - i) = -E.

,

i=

1,2,...,n

(63)

-= B n ( i ) + D n B n ( n + 1 - r - i) = On

,

i=l,2,...,n-r

= -Bn(i) + E n B n ( n + l - r - i) = E,

,

i = l,2,... ,n - r

(64)

and

{

~,+a(i) B n + l ( n + 1 - r)

(65)

In spite of that the description of the proposed algorithm seems complicate, calculations based on it are relatively simple. In fact, first for n' = 1 , 2 , . . . , r matrices A ~ , ( i ) are determined by the Levinson-Durbin algorithm, after that, matrices A t ( i ) will be calculated as below: -

determine ~x(1) by (47), and Ha by (50);

using the recursive formula (51) calculate ~2(i),... and so on we get At(i) - ~r(i). Furthermore for n = r + 1,r + 2 , . . . , N, first we calculate Dr and Er by (58) and (61), respectively, then using the recursive formulae (62)-(65) we compute A r + l ( i ) , A r + , ( i ) , Br+a(1) and Br+l(1). After that, it is neceessary to come back to calculate Dr+a and Er+a before repeating the computations of the (r + 2)-th order matrices,.., and so on we continue such recursive calculations upto the desired order N. -

4. THE RELATION BETWEEN THE PARAMETERS OF THE AR AND ARMA MODELS The AR model is a special case of the ARMA model, when the order p of the MA part in Eq. (2) is equal to zero, or in accordance with it the operator function B ( z ) in the definition (8) is the unit matrix: B ( z ) = I.

Recursive estimation of ARMA parameters

299

In this section we will show how the solution of the Yule-Walker equation (18) becomes the one for the AR model, i f p = 0. In this case Eq. (18) will be written as NAR

A(z)C(k)-

I-

E

(66)

M i z - i ~ C(k) = 0 .

i=1

/

We can see that this equation still differs from the Yule-Walker equation of the AR model by the term

NAR

i~

.._~ M i z - I C(k). So a question arises that this term vanishes or not,

I-

/

and if yes, then under what conditions. In Kishida (1982, 1984) the innovation (3) was used during the construction of the ARMA model from the Langevin physical equation. This innovation was rewritten by us in section 2 as the expression (11). From (4) and (16) the term M Y ( k - 1) now can be expressed by y(k-1)

MY(k

1)

.

.

[M,,M2, .

.

.,MNAR] .

}

y ( k - 2)

= NAa F , M,y(k - i )

.

(67)

i=l

y(k -- NAR) So the innovation (3) becomes: NA R

(68)

v(k) = u(k) - y ~ M,u(k - i) i=1

To find the meaning of the innovation (3) and (68), we have to clarify the meaning of Mi. This is done below. T h e o r e m : The matrices Mi given in (16) are the solution of the Yule-Walker equation of the AR model with the NAn-th order. P r o o f i For the theorem, we will prove that

I

NAR

I-- E

Miz-i

1

C(k)=O ,

that is

i=1

NA~

MiC(k-i)=C(k)

,

(69)

i=1

k = 1, 2 , . . . , NAR . In fact, let us introduce a matrix ~ as

=

c(0) C(-1) C(-2)

\c(-g

iA

R-t-l)

c(1) C(O) C(-1)

... ... ...

" " ' " C(--NAR+2) ...

C ( N A R - 1) C(NAR-- 2) C(NAR-- 3) " C(O)

= [~I,~2,'-',~NAR],

(70)

300

TRAN DINH TRI

where ~bj is the j - t h block column of q:

CJ =

C(j .- 2) C(j -" NAR)

We can see that q is a symmetric matrix, and its block transposition yields a so called block Toeplitz matrix• The matrix M in Eq. (16) now can be rewritten as M = [C(1),C(2),...,C(NAR)] kv-1 , then MqJ = [C(1), C ( 2 ) , . . . , C(NAR)]

(71)

On the other hand, using the definition (70) of qJ we have: M ~ --= [M¢I , M ¢ 2 , . . . , M ¢ N a n ] •

(72)

From (71) and (72) we can see that [MO1,M¢2,...,MCNAR]-- [C(1),C(2),...,C(NAR)]

.

(73)

Using the expression (17) we have

M C j ---- [M1, M2, • .., MNA•]

C(j •- 2)

= NAR Z M i C ( j - i) .

(74)

i=1

C(j - NAR) Substituting (74) into (73) we obtain the result:

~ --

MiC(1-

i), Z

MiC(2-

i),..., Z

i=l

MiC(NAR-

i)

= [C(1),C(2),...,C(NAR)].

i=1

(75) This result allows us to conclude that NA R

M~C(k

- i) :

C(k)

,

k = 1,2,...,

NA,

•

(76)

i=1

That is what we have to prove. Now based on the theorem proved above we can show that the number NAR is the order of the AR model (in accordance with this the subindex A R is attached to the nmnber NAR), then it can be chosen by the AIC (Akaike's Information Criterion), the F P E (Final Predition Error) or the MEM (Maximum Entropy Method) (Kishida, 1980). Furthermoro,

Recursiveestimationof ARMAparameters

301

for n < NAR the Yule-Walker equation (18) for the ARMA model in the case p = 0 exactly agrees with the one for the AR model, and of course in this case the ARMA solution, which is obtained by the algorithm proposed in the previous section, becomes the AR solution obtained by the Levinson-Durbin algorithm. Finally, with the help of the above theorem we have also shown an important meaning of the innovation ~/(k) defined by the expressions (3) or (68) that this is not else than the residual noise of the AR model. 5. CONCLUDING REMARKS On the basis of the algorithm proposed in section 4, we have built up a progrmn package called ARMAANA to evaluate the P W R and BWR reactor noise diagnostic measurement data. The application of this program package for different pratical data sets is now being continued with the aim to evaluate how better the results of the ARMA model are in comparison with the ones of the AR model. The first results obtained in these tests are encounaging and will be presented in forthcoming communication.

REFERENCES

Fukunishi K. (1977) NucI. Sci. and Eng., 62, pp. 215-225. Kishlda K. and Sasakawa H. (1980) J. Nucl. Sci. Technol., 17, 16. Kishida K. (1982) Physical Review A, Vol. 25, No. 1, pp. 496-507. Kishida K. (1984) Y. Math. Phys., Vol. 25, No. 5, pp. 1308-1313. Kishida S. K. and Sekiya T. (1976) J. Nulc. Sci. Technol., 13, 19. Tomita K. and Tomita H. (1974) Prog. Theor. Phys., 51, 1731. Tran Dinh Tri (1990) KFKI Report, 1990-45/(], Budapest, Hungary. Upadhyaya B.R. et al. (1980) Ann. nucl. Energy, 7, pp. 1-11. Whittle P. (1963) Biometrika, 50, 129.