HUNGARIAN JOURNAL OF INDUSTRIAL CHEMISTRY VESZPREM Vol. 30. pp. 215- 218 (2002) TIME SERIES AND THE ALGEBRAIC MATRIX RICCATI EQUATION C. STOREY (The Institute of Mathematics and Simulation Sciences, Faculty of Computing Sciences and Engineering, De Montfort University, The Gateway, Leicester, LEl 9BH, UK) Received: September I 0, 2002 One way of modelling certain kinds of time series is via the Yule-Walker equations. These are a set of (over-determined) linear equations for estimating the parameters in the models. The coefficients in these equations are estimates, obtained from the data, of the autocorrelations. In this paper two ways of solving the Yule-Walker equations are considered. The first is the well known method using the pseudo-inverse and the second uses the algebraic matrix Riccati equation. A number of numerical examples are used to illustrate and compare the two different approaches. Keywords: autoregressive time series, Yule-Walker equations, algebraic matrix Riccati equations, total least squares Introduction Determination of the autoregressive parameters ai,i = O,l,···,n, (1) in the stochastic process (2) where (3) i=o is investigated. In Eq.(2), y(t) denotes noisy data, w(t) is white noise and D(q-1 ) is a polynomial containing the moving average parameters of the process. The expression q-1 in Eq.(2) is the unit delay operator defined by q -1 y(t) = y(t - 1). More extensive details and background material can be found in Ref. [1]. Here the Yule-Walker method is used to estimate the parameters ai by two different techniques. Firstly the direct method involving the inverse or pseudo-inverse of the autocorrelation matrix, depending on the invertability or otherwise of the latter, is used. Secondly an algorithm is used which is based on the result that a total least squares solution to the Yule- Walker equations can be found by solving an appropriate matrix Riccati equation. A number of numerical examples from Ref. [1] serve to illustrate the two techniques and to compare the results with a more sophisticated stochastic analysis given in Ref. [1]. A simil'!r investigation has been carried out in Ref. [5] but the TSL solutions to the Yule~ Walker equation are obtained by singular value decomposition and simulated time series are used. The Yule-Walker Method For a pure (i.e., D(q-1 )= 1) AR process the Yule- Walker equations are (Ref.[ 1 ], Ref. [2]), r rli Tn li ... rn-I a I 'i To ···~z-2 a'2 rJ = (5) rm rm-1 ••• rm-(n-1) an r., 216 Table 1 Numerical results for Series C nxm lxl lxlO lx50 lxlOO TENT FINAL LS -.805 -.814 -.847 -.881 -.81 -.82 TLS -.805 -.815 -.854 -.892 -.81 -.82 Table 2 Numerical results for Series D nxm lxl lxlO lx50 lxlOO TENT FINAL LS -.861 -.860 -.908 -.908 -.86 -.87 TLS -.861 -.860 -.916 -.912 -.86 -.87 ck is an estimate of the relevant unknown autocovariance given by ck =_!_ ~(Yt-Y y Yt+k- y), k=O,l,2,··· N t=! A and N is the number of observations in the time series. In case m > n then the linear equations (Eq.(S)) are overdetermined and the pseudo-inverse has to be used. Numerical Investigation Four different autoregressive time series taken from Ref. [1] have been used as numerical illustrations. ~ample I Series C This is a time series, with N = 226 observations, of temperatures in a chemical process taken at one minute intervals. A model with two possible sets of parameters, one tentative (tent) the other more finalized (final) is given in Ref. (1). These are the following: Tent: Vy(t)- 0.81Vy(t- I) = w(t) , Final: Vy(t)-0.82Vy(t-l)=w(t). The results of the numerical investigation are given in Table 1. The table entries are solutions (rounded to 3 decimal places) of n x m Yule~ Walker equations by the least squares method (LS) and by the total least squares method (TLS). The former is found using the pseudo- inverse matrix and the latter using Algorithm l of [3]. (A convenient check on the LS solution is also obtained by using Al as explained in [4]). For this example a closer analysis shows that m = 22 is about the best number of equations to use and TLS is slightly better than LS. Example 2 Series D The tentative and final models for this example are: Tent: y(t)- 0.86y(t- l) = '"1t) Table 3 Numerical results for Series E. Model I nxm 2x2 2x10 2x50 2x99 TENT FINAL LS -1.318 -1.398 -1.475 -1.515 -1.32 -1.42 .634 .704 .724 .770 .63 .73 TLS -1.318 -1.432 -1.549 -1.586 -1.32 -1.42 .634 .733 .790 .834 .63 .73 Table 4 Numerical results for Series E. Model2 nxm 3x3 3x10 3x50 3x99 TENT FINAL LS -1.369 -1.821 -2.104 -2.155 -1.37 -1.57 0.740 1.378 1.775 1.849 0.74 1.02 -0.0805 -0.367 -0.571 -0.599 -0.08 -0.21 TLS -1.369 -1.975 -2.363 -2.412 -1.37 -1.57 0.740 1.609 2.188 2.277 0.74 1.02 -0.0805 -0.485 -0.783 -0.825 -0.08 -0.21 Final: y(t)-0.87y(t-l)=w(t). There are N = 310 observations and the data is a set of chemical process viscosity readings taken at hourly intervals. The numerical results are shown in Table 2. For this example m = 19 appears to be the best number of equations to take with TLS getting nearer to - 0.87 than LS with the same number of equations. Example 3 Series E For this, much studied example, the data shows the number of sunspots that occurred in each year from 1770 to 1869; so there are 100 observations. Two different sets of parameter values are given for each of two different suggested models: Modell Tent: y(t) -1.32y(t -1) + 0.63y(t- 2) = w(t) , Final: y(t) -1.42y(t -1) + 0.73y(t- 2) = w(t) , Model2 y(t) -1.37y(t -1) + Tent: , +0.72y(t-2)-0.08y(t-3) = w(t) Final: y(t) -1.57 y(t -1) + . + 1.02y(t- 2)- 0.21Y(t- 3) = w(t) T abies 3 and 4 give the numerical results in this case. For the first model m = 9 or 10 give good approximations to the parameter values given in [ 1] and for the second model m = 7 or 8 seem to be reasonable values. Table 5 Numerical results for Series F nxm 2x2 2x10 2x50 2x69 TENT FINAL LS 0.32 0.319 0.331 0.365 0.32 0.34 -0.18 -0.170 -0.178 -0.149 -0.18 -0.19 TLS 0.32 0.323 0.385 0.451 0.32 0.34 -0.18 -0.168 -0.157 -0.104 -0.18 -0.19 Example 4 Series F Here the time series F consists of 70 readings of the yields from a batch chemical process and the suggested model is: Tent y(t) + 0.32y(t-1) -0.18y(t- 2) = w(t) Final: y(t) + 0.34y(t -1)- 0.19y(T- 2) = w(t) The results of the calculations are shown in Table 5. For this example the second parameter increases with m so the best m for this parameter is 2! For the first parameter, however, m = 25 gives a value of -0.338 with 0.169 for the second parameter. "Portmanteau" lack-of-fit test A test for model adequacy, which depends on the first K auto-correlations of the residuals rk ( Q) , is proposed in [1]. The statistic K Q = n(n+2) L.rk2(t1)/(n-k) (6) k==l is computed and referred to a table of the chi-squared distribution. In Eq.(6) n = N -d where d is the degree of differencing required to approximate to stationarity in the time series being analysed. For series D, with K = 25 , there are 24 degrees of freedom and Eq.(6) becomes 25 Q =310x31iL,rk\t1)/(310-k) , k=I which is readily computed by MATHEMATICA, for example, as follows: Subtract the mean of series D (9.13258) from each of its elements and can the resulting series, augmented by a back-forecasted value as its initial value, s(D). Compute the series of residuals. s(D) =Table [Part [s(D),z1- $ Part[s(D),i-l],{i,2,310}] with $ the estimate of the model parameter. Compute Q as follows: s(D)=Table[Part[s(D),z1-0.0024394l,{i,l,310}] (0.00243941 is the mean of the previous s(D)). ro(a) = Sum[Part[s(D),i]* Part[s(D),i],{i,l,310}] 217 g = Table[Sum[Part[s(D),i] * * Part[s(D),i+ j],{i,l,310- j}],{j,l,l}]lr,,(a) (This produces a table of l auto-correlations of residuals) Q = 310*312Sum[Part[g,if2/(310-i),{i,l,25}]. Taking $ = -0.87 (the final value from [1]) gives Q=l1.5 Now when m x n = 19xl the parameter obtained by LS is 0.868034 and by TLS it is 0.868979. The corresponding values of Q are 11.44 and 1 1.47. Similarly when m x n =20xl the values for are 0.870544 and 0.87177 with corresponding values of 11.51 and 11.55 for Q. Thus TLS gets closer to the value of Q given in [1] with fewer linear equations. Conclusion A numerical study has shown that the TLS method is slightly superior to the LS method for finding the unknown parameters in models of auto~regressive time series. The techniques used can be extended to estimation of moving average parameters [2]. Special purpose methods for solving the Yule-Walker equation are also discussed in [2]. Finally, when Algorithm 1 is used on its own, and not in comparison with other methods as it is here. the computations should be stopped as soon as some test of lack of fit (such as the portmanteau test from [1]) is adequately satisfied. ai A,D q·l y{t) w(t) n m y N d SYMBOLS autoregressive parameters polynomials in q-1 unit delay operator noisy data white noise number of autoregressive parameters number of Yule-Walker equations mean ofy(t) number of observations in time series degree of differencing for stationarity Q statistic used in Portmanteau test i, j, k running incides l number of autocorrelations taken in Portmanteau test. REFERENCES 1. Box G. E. P .• JENKINS G. M. and REINSEL G. C.: Time Series Analysis: Forecasting and Control. 218 Prentice Hall Inc., Englewood Cliffs, New Jersey, 3rd Edition, 1994 2. SODERSTROM T. and STOICA P.: System Identification, Prentice Hall International, 1989. 3. STOREY C.: Hung J. Ind. Chem.J997, 25,305-308 4. STOREY C.: Ibid, 2001,29,67-70 5. STOIAN A., STOICA P. and VAN HUFFEL S.: Scientific Bulletin, Electrical Eng. Series, Polytechnic Institute of Bucharest, 52, No.2, 81-89 Page 217 Page 218 Page 219 Page 220