- Email: [email protected]
CONTINUOUS-TIME AUTOREGRESSIVE SPECTRAL ANALYSIS FOR IRREGULARLY SAMPLED DATA
14th IFAC Symposium on System Identification, Newcastle, Australia, 2006
CONTINUOUS-TIME AUTOREGRESSIVE SPECTRAL ANALYSIS FOR IRREGULARLY SAMPLED DATA Piet M.T. Broersen Department of Multi Scale Physics, Delft University of Technology
Abstract: A continuous-time autoregressive spectral estimator is introduced that applies the principles of a discrete-time automatic equidistant missing data algorithm to unevenly spaced data. This time series estimator approximates the irregular data by a number of equidistantly resampled missing data sets, with a special nearest neighbor method. The ARMAsel-irreg algorithm estimates and automatically selects a discrete-time AR model from a number of candidates. This selected model often has a number of spurious high frequency poles, which are incompatible with the continuous character of the irregularly sampled signal. Those spurious poles can be eliminated, by transforming only the poles of the discrete time model with a positive real part to matching continuous-time poles. The estimated continuous-time spectra can be accurate at frequencies much higher than the mean data rate. Copyright © 2006 IFAC. Keywords: autoregressive model, nearest neighbor resampling, slotting, spectral estimation, time series analysis, uneven sampling, order selection, parametric model.
1. INTRODUCTION
normalization reduces the variance of the estimated autocovariance function. Fuzzy slotting produces a smoother autocovariance by distributing products over multiple time slots. Autocovariances estimated by slotting techniques are not guaranteed to be positive semi-definite. This results in spectra that can become negative at frequencies where the power is weak. Slotting methods require very large data sets.
Meteorological data, astronomical data or turbulence data obtained by Laser-Doppler anemometry are often irregularly sampled, due to the nature of the observation system. This has the advantage that the highest frequency that can be estimated is higher than half the mean data rate, which is the upper limit for equidistant observations. Many estimation techniques exist for unevenly spaced data. Lomb (1976) and Scargle (1989) described a Fourier estimator. Bos, et al. (2002) compare several variations of slotting or resampling methods. Larsson and Söderström (2002) replace the derivative operator in the continuous-time model by a discrete time approximation and apply it to a second order process. Lahalle, et al. (2004) estimate continuous time ARMA models, requiring the explicit use of a model for the irregular sampling instants. Larsson and Larsson (2002) give the Cramer-Rao lower bound for continuous-time ARMA models.
Resampling techniques reconstruct a signal at equal time intervals. After resampling, the equidistant data can be analyzed using the periodogram or time series models. Spectral estimates at higher frequencies will be severely biased. Adrian and Yao (1987) described Sample and Hold reconstruction as low-pass filtering followed by adding noise. These effects can in theory be eliminated using the refined Sample and Hold estimator (see Benedict, et al., 2000). In practice, all spectral details smaller than the bias are lost. Nearest Neighbor resampling has similar characteristics. The spectra are strongly biased for frequencies higher than about 20 % of the mean data rate. The noise and filtering effects of equidistant resampling set limits to the achievable accuracy of resampling methods. This precludes the accurate estimation of spectra at higher frequencies where the resampling noise blurs details,
Benedict, et al. (2000) give an extensive survey of methods for irregular data. Slotting methods estimate an equidistant autocovariance function from irregular samples. Slotting has been refined with normalization and fuzzy slotting. Local 849
smaller than the bias and hides spectral slopes.
any slot. Therefore, multi shift resampling is used, where M different equidistant missing data signals are extracted from one irregular data set and all small slots of width w are connected in time.
Bos, et al. (2002) introduced an estimator that can be perceived as searching for uninterrupted sequences of data that are almost equidistant. The contiguous sequences of different lengths can be analyzed with an irregular version of the Burg (1967) algorithm for segments. Slotted nearest neighbor resampling creates an equidistantly sampled signal by replacing an irregular observation time instant by the nearest equidistant resampling grid point, but only if it is within half the slot width. The slot width is equal to the resampling time Tr. Some grid points may remain empty. A smaller slot reduces the bias. The bias of slotted nearest neighbor is very much smaller than the bias without slotting. The reason is that one original irregular observation can never appear at multiple resampled time instants if slotting is applied. A disadvantage of this slotted Burg resampling method is that still very large data sets are required to obtain some uninterrupted equidistant sequences of sufficient length for the irregular Burg algorithm. It turned out that a non-linear maximum likelihood algorithm for missing data, developed by Jones (1980), could give a better solution, also if a limited amount of data is available.
This paper studies the possibility of automatic AR order selection for irregular data. It also describes how the transformation of the discrete-time model, with the automatically selected order, to a continuous-time model is favorable for the spectral accuracy. By simply eliminating poles of the selected estimated discrete-time model that cannot possibly be produced by poles of the continuous-time data, an accurate autoregressive spectrum is obtained. 2. TIME SERIES MODELS A discrete-time AR(p) model can be written as (Priestley, 1981) xn a1 xn 1 a p xn p
(1)
Hn ,
where Hn is a purely random process of independent identically distributed stochastic variables with zero mean and variance VH2. Assume that data represent a stationary stochastic process. The power spectral density h(Z ) of an AR(p) model is completely determined by the parameters in (1) together with the variance VH2 and is given by:
To study the possibilities of discrete-time models for irregular data, a separate study has first been made to the simpler problem of equidistant observations with missing data. In that case, Jones (1980) described an efficient method to calculate the true likelihood. A survey of existing missing data methods and a numerically stable version of the maximum likelihood algorithm for autoregressive models of missing data problems has been given by Broersen, et al. (2004a, 2004b). The performance of that automatic time series algorithm outperforms all other methods for missing data problems. The best method for missing data will be applied to irregular data.
h(Z )
V H2 1 2S A e jZ i p
2
V H2 2S
1 p
2
.
(2)
1 ¦ ai e jZ i i 1
Furthermore, the invariance property of maximum likelihood theory shows that the transformation (2) of efficiently estimated time series parameters also gives an efficient estimate of the power spectrum. Broersen (2002) developed a computer program to estimate the parameters and to select the model order with no other input than the data. It requires no user interaction to determine the best model for a given set of random observations.
Broersen and Bos (2006) showed that the continuous time maximum likelihood approach of Jones (1981) mostly fails to converge for irregular observations. If more than 2 parameters have to be estimated, the numerical surface of the continuous likelihood is too rough for global convergence, due to local maxima. Modifications are required to apply the automatic algorithm for equidistant missing data of Broersen, et al. (2004a,b) to irregularly sampled data. That needs multi shift slotted nearest neighbor resampling and uses a slot width that is only a fraction of the resampling time. This produces several equidistant data sets, with slightly different starting points that are shifted over the slot width. Time instants are fixed to an equidistant resampling grid where the original irregular sampling instant is never further away from the grid point than half a slot width. The resampling grid time Tr determines the discrete-time frequency range. By taking a slot width w that is only a fraction 1/M of the resampling time Tr, the bias is further reduced. This, however, gives disjoint intervals where many irregular times ti are not within
A similar automatic, approximate maximum likelihood, ARMAsel-mis program has been developed for the equidistant missing data problem. The accuracy of the ARMAsel-mis spectra is better than the spectra obtained with many other methods from the literature, in simulations with missing data. Broersen, et al. (2004b) have given examples where the estimation of time series models in missing data problems was efficient, meaning that the accuracy of the resulting model approached the limit of the achievable accuracy. The application of that discretetime missing data program to continuous-time irregular data will be studied here. That is necessary because the irregular maximum likelihood method of Jones (1980) mostly fails to converge and no other reliable automatic spectral estimation methods for finite or small samples of irregular data are available. The missing data AR order selection criterion GIC ( p )
850
LH D p,
(3)
uses the minimum negative log likelihood LH with a penalty factor D that depends on the missing fraction, 3 for less than 25 % missing, 5 for more than 75 % missing and 4 in the range between. T0 /MTr approximately gives the remaining fraction where T0 equals 1/f0, which in turn is the mean data rate.
the integral of the squared difference between the logarithms of spectra can be computed for an arbitrary frequency range 1/2Ta :
A continuous-time AR(p) model can be written as (Priestley, 1981) (4) s p D1s p1 D p x(t ) H (t ) , sx(t ) dx dt where s denotes the differentiation operator and H(t) is white noise with as correlation a G function and as spectrum the constant V 2. The power spectral density h(Z ) of an AR(p) model is given by the parameters :
The hat indicates a spectral estimate. By limiting the integration to the definition area 1/2Tr of the discrete spectrum, it is possible to attribute a number to the accuracy of discrete-time or continuous-time approximations.
h(Z )
V2 2
p
¦D
i
e
, D0
1, f Z f .
SD
NTa 2S
S / Ta
³ S
/ Ta
2
^
`
ªln ^h(Z )` ln hˆ(Z ) º dZ . ¬ ¼
(6)
3. ARMASEL FOR IRREGULAR DATA In the first application of the ARMAsel-irreg algorithm, Broersen and Bos (2006) reported good results for known low order models. Automatic order selection sometimes failed because of the peculiar behavior of the discrete-time likelihood function as a function of the model order. Often, spurious high frequency peaks seem statistically significant. That is illustrated in Fig. 1 for an example with AR(1) data.
(5)
jZ ( p i )
i 0
The discrete-time spectrum in (2) covers the frequency range Z between -S and S; the continuoustime spectrum has Z between -f and f. Transformations between models must always be a compromise of properties. Franklin, et al. (1998) describe transformations from a control point of view. The bilinear or Tustin transformation is often a good choice for the conservation of filter properties. A disadvantage is that the frequency scale is warped. The zero order hold equivalents can give the same outputs for continuous-time and discrete time systems at the sampling instants. Söderström (1991) gives several solutions for the computation of stochastic continuous-time models from discretetime ARMA models. He showed that negative real discrete-time poles and also some more specified areas of zero positions do not give any equivalent continuous-time solution. Those stationary and invertible ARMA models cannot possibly be obtained by sampling a continuous-time process.
Fig.1. Spectrum of discrete-time AR models, shifted downwards with a factor 10 for each new order. 8 AR models are estimated from N=100 irregular samples of the true continuous AR(1) process with the pole at - 0.1S, resampled with Tr = 1/2f0, where f0 =1 is the mean data rate; w = Tr. In this example, the AR(3) model was selected.
For spectral analysis, the conservation of the location of spectral peaks and valleys is the most important. That can be obtained with a technique called pole and zero matching. A pole at s = D in (4) will give a pole at exp(DTr) in (1). A pole at s = D + iE will give a pole at Uexp(iM) in (1), with U = exp(DTr) and M = ETr, (Franklin, et al., 1998). Stationary or stable processes must have only poles with a negative real part in the continuous s-domain or a radius less than 1 in discrete time models. The same matching rules can be used for the inverse transformation of discrete-time to continuous-time models. All discrete poles within the unit circle have a stable transform, except negative real poles for which no sensible matched pole transform can be defined. AR processes have no zeros, otherwise the same rules would be used for the matched transform of zeros.
The order 3 was selected in Fig. 1, although the best approximation was given by the AR(1) model. It is curious that the AR(3) model has a strong peak at the end of the frequency interval. This belongs to a negative real pole, which cannot belong to a sampled continuous-time process. All the higher order models show spurious details. Strange peaks can mostly be found in the range between 1/4 Tr to 1/2Tr. What is worse, models with such ‘impossible’ peaks at high frequencies are often selected by GIC(p) of (3).
Accuracy measures have been defined by de Waele and Broersen (2000) that can compare the quality of the discrete-time model with true or with aliased continuous spectra. With the spectral distortion SD,
Simulations with higher order generating processes have a similar behavior. Estimated low order models may have acceptable spectra. But at and above the true order, and often already below, estimated
851
models often have spurious peaks that seem to be statistically significant. Most peaks appear at the highest discrete frequency 1/2Tr, many in the last half of the frequency range, above 1/4Tr, and only very few are found for lower frequencies. Order selection criteria often select one of the discrete-time models with a peak at 1/2Tr, a model without a continuous time equivalent. It has been verified that this agrees with the minimum of the log likelihood in those examples. In other words, the likelihood decreases with a statistically very significant amount by adding a negative real pole, giving a spectral peak at the frequency 1/2Tr. No value of the penalty factor D in (3) would give a good performance of the order selection criterion. An analysis of the bias of ARMAsel-irreg by Broersen and Bos (2006) showed that the frequency range between 1/4 Tr to 1/2Tr suffers from serious bias in examples with strong spectral slopes at those frequencies. That frequency range is very sensitive to bias and to spurious peaks.
4. PERFORMANCE OF THE NEW ESTIMATOR Simulations with a known spectrum are a first step in testing new algorithms. Test data was generated with N irregular continuous-time data points of a known AR process. The resulting irregular data has time intervals that are Poisson distributed. In simulations, the true properties of the data are known. Hence, the quality of estimated results can be established with the SD of (6). The vertical scales in the spectral plots are sometimes shifted, to increase the visibility.
An example can be given where a missing data problem with a pole at a1 could have a minimum of the likelihood for a pole at – a1. That can happen for an AR(1) process with parameter a1 if the data are missing with a certain pattern. Suppose that all odd observations are missing, only even time instants remain, for T = 1. The autocorrelation function of this process is the series (-a1)k and that is identical for the parameters + a1 or – a1, for all even values of k. The convergence of the optimization algorithm for the AR(1) parameter depends on the starting values then, because the likelihood is symmetric. That has been verified. The same regular pattern of only even samples could theoretically also be found for irregular sampling. It gives at least an explanation that it is possible to find those strange, ‘impossible’, peaks in an algorithm without programming errors.
Fig. 2. True spectrum, discrete-time ARMAsel-irreg and two continuous-time estimates for a true AR(2) spectrum. Tr = 1/4f0, with f0 = 1 and slot width w = Tr /2. N = 250. ARMAsel-irreg selected AR(4). Fig. 2 gives spectra for 250 irregular observations of an AR(2) process with 1 peak. The ARMAsel-irreg program selected the discrete-time AR(4) model from candidates AR models until order 6, with a log likelihood value LH equal to - 152.8. The estimated discrete-time AR(2) model had LH= - 122.3. That spectrum fitted much better than the selected AR(4) model. However, after matched transformation to continuous-time, the ARMAsel continuous Matched AR(4) model had a better spectral estimate than the transformed AR(2) model. The matched continuous model eliminates all discrete poles with negative real parts, which agrees with frequencies greater than Tr/4. The bilinear continuous transformed spectrum keeps the wrong peak of the ARMAsel-irreg AR(4) model and shifts it to a new frequency in Fig. 2. This shifting property as well as the conservation of the peak with a frequency above Tr /4 makes the bilinear transformation less adequate for this application. Repeating this simulation 10 times always gave a selected AR(3) or AR(4) model, with about the same accuracy for the final ARMAsel continuous matched model. The true continuous process has a spectrum with an infinite frequency range, the estimated discrete-time model has a frequency range until 2f0 in Fig. 2. The range of the transformed continuous-time models is also infinite. However, as those models are derived from the discrete-time model, it is an open question how reliable high frequencies are. An
Two completely different algorithms for the log likelihood for missing data have been programmed independently by Broersen, et al. (2004a, 2004b): the exact likelihood as formulated by Jones (1980) and an approximation. Both algorithms lead to the same estimates for the likelihood and to the same spurious high frequency spectral peaks. Broersen and Bos (2006) also showed that the continuous likelihood has an undesirable character, a conclusion that was also supported by other authors. It still is an open question whether the regular even sampling pattern of the previous paragraph is a possible explanation for peculiar behavior in Fig. 1. No other explanation is available yet. A practical solution for the order selection problem for irregular observations will be investigated by switching to continuous time models. First, ARMAsel-irreg selects the discrete-time model with the minimum of GIC(p) of (3). Afterwards, all poles of the selected model with frequencies above 1/4Tr will be eliminated before a matched transform is applied. By choosing the value of Tr small enough, it can be made sure that no true spectral peaks are in that frequency range. In this way, the frequency range that is sensitive to bias and to spurious spectral peaks is eliminated in the continuous-time model.
852
Fig. 4 gives the result for an AR(4) process with two strong spectral peaks. All discrete-time higher order models showed a strong peak at the end of their frequency range. In this example the AR(4) model would also give a good spectrum. In another run with 10000 observations, the estimated AR(4) model found only the second true peak, and an extra peak at the end of the frequency range. It missed the first true peak. Obviously, the spurious high frequency peak was more significant in the likelihood. In that run, the selected model order was much better than the true order. That often has been seen in individual simulation runs. Knowing the true order mostly has no advantage for irregular data in our simulations. Using selected orders is generally more accurate.
intermediate range is used for the visual representation, to indicate the differences. The perfect match of the estimated Matched spectrum at the high frequencies is caused by the agreement of the number of true and estimated poles. Each pole gives a high frequency spectral asymptote of –20 dB per frequency decade in continuous modeling.
It has been tried to use the parameters of the selected continuous matched ARMAsel model as starting values for a further search of a continuous maximum likelihood model. The very irregular surface of the continuous likelihood function, as reported by Broersen and Bos (2006) prevented this approach. Sometimes, the estimated model did not improve but remained the same. It could also happen that the continuous likelihood converged to a much poorer model than the starting position. It never occurred in the simulations that the AR(4) example converged to a much better model. This continuous maximum likelihood post processing could only be used fruitfully for AR(1) and AR(2) examples, where also the approach of this paper did not give any problem.
Fig. 3. True spectrum, ARMAsel-irreg and two identical continuous matched estimates for AR(3) data. Tr = 1/2f0 with f0 = 1 and slot width w = Tr /2 , N = 1000. ARMAsel-irreg selected AR(3). Fig. 3 gives the result for a smooth AR(3) process with the three real poles at -0.1S, -0.4S and -S. AR(3) was selected with ARMAsel-irreg and higher order models showed only a moderate tendency for spurious details. In general, examples with strong details in the spectra have a tendency to overestimate the order and to have more spurious peaks for the same sample size. In this example, the selected and the true order model coincided because AR(3) was selected. Remind the vertical shift of other spectra.
Fig. 5. True, discrete-time ARMAsel-irreg spectrum, and a continuous matched estimate for a true AR(5) spectrum. Tr = 1/2f0 with f0 = 1 and slot width w = Tr /2. N=2000. ARMAsel-irreg selected AR(3). Fig. 5 has a rather smooth true spectrum with 5 real poles at -0.1S, -0.2S, -0.4S, -0.8S, and -1.6S. AR(3) was selected by ARMAsel-irreg. The estimated AR(5) model had one negative real pole, which was eliminated. The different orders can be recognized in the different asymptotic slopes, -60, -80 and –100 dB for AR(3), AR(4) and AR(5), respectively. The performance of this estimator and the spurious peaks for greater sample sizes still have to be investigated.
Fig. 4. True, discrete-time ARMAsel-irreg spectrum, and a continuous matched estimate for a true AR(4) process. Tr = 1/4f0 with f0 = 1 and slot width w = Tr/2. N = 1000. ARMAsel-irreg selected AR(6).
853
Autoregressive spectral estimation by application of the Burg algorithm to irregularly sampled data. IEEE Trans. on Instrumentation and Measurement, 51, pp. 1289-1294. Broersen, P.M.T. (2002). Automatic spectral analysis with time series models. IEEE Trans. on Instrument. and Measurement, 51, pp. 211-216. ARMASA Matlab toolbox, freely available at www.mathworks.com/matlabcentral/fileexchange.
Broersen, P.M.T. and R. Bos (2006). Estimation of time series models from irregularly spaced data. IEEE Trans. on Instrumentation and Measurement, offered for publication. Broersen, P.M.T., S. de Waele and R. Bos (2004a). Application of autoregressive spectral analysis to missing data problems. IEEE Trans. on Instrument. and Measurement, 53, pp. 981-986. Broersen, P.M.T., S. de Waele and R. Bos (2004b). Autoregressive spectral analysis when observations are missing. Automatica, 40, pp. 1495-1504. Burg, J. P. (1967). Maximum entropy spectral analysis. Proc 37 th Meeting Soc. Of Exploration Geophysicists, 6 pp. Franklin, G.F., J.D. Powell and M. Workman (1998). Digital Control of Dynamic Systems. Addison Wesley, Menlo Park, CA, USA. Jones, R.H. (1980). Maximum likelihood fitting of ARMA models to time series with missing observations. Technometrics, 22, pp. 389-395. Jones, R.H. (1981). Fitting a continuous time autoregression to discrete data. Applied Time Series Analysis II, (Ed. D.F. Findley), pp. 651682. Lahalle, E., G. Fleury and A. Rivoira (2004). Continuous ARMA spectral estimation from irregularly sampled observations. Proc. IEEE/IMTC Conf., Como, pp. 923-927. Larsson E.K. and T. Söderström (2002). Identification of continuous-time AR processes from unevenly sampled data. Automatica, 38, pp. 709-718. Larsson, E.K. and E.G. Larsson (2002). The CRB for parameter estimation in irregularly sampled continuous-time ARMA systems. IEEE Signal Processing Letters, 11, pp. 197-200. Lomb, N.R. (1976). Least squares frequency analysis of unequally spaced data. Astrophysics and Space Science, 39, pp. 447-462. Priestley, M. B. (1981). Spectral Analysis and Time Series. Academic Press, London, U.K. Scargle J.D. (1989). Studies in astronomical time series analysis III. Fourier transforms, autocorrelation functions, and cross-correlation functions of unevenly spaced data. The Astrophysics Journal, 343, pp. 874-887. Söderström, T. (1991). Computing stochastic continuous-time models from ARMA models. Int. J. Control, 53, pp.1311-1326. Waele S. de and P. M. T. Broersen (2000). Error measures for resampled irregular data. IEEE Trans. on Instrumentation and Measurement, 49, pp. 216-222.
Fig. 6. True, discrete-time ARMAsel-irreg spectrum, and two continuous estimates for a true AR(6) spectrum. Tr = 1/4f0 with f0 = 1 and slot width w = Tr /2 , N = 2000. ARMAsel-irreg selected AR(10). The example in Fig. 6 demonstrates that sometimes the spurious peaks are estimated in AR(6) before the true peaks. The elimination of the discrete-time poles with negative real parts from AR(10) is very favorable here. AR(10) was the highest candidate.
5. CONCLUSIONS A new estimator is introduced that first fits discretetime AR models to multi shift slotted nearest neighbor resampled segments from irregular data. The ARMAsel-irreg algorithm automatically selects the AR order. Afterwards, the reliable poles of the selected discrete-time AR model are transformed to matched continuous-time poles. Only estimated poles with a positive real part are used. All poles with negative real parts are eliminated, because they mostly belong to spurious spectral details. In simulations with few irregular data, the results are much better than what can be obtained with any other known existing spectral estimation technique. Multi shift slotted nearest neighbor resampling will give accurate spectra, with a small bias, if the dynamic spectral range is limited. For a large dynamic range, a small slot width reduces the expected bias. The elimination of continuous-time poles with a frequency above 1/4Tr further reduces the bias. In practice, Tr can be chosen small enough to make sure that no true process poles will be eliminated. REFERENCES Adrian R. J. and C. S. Yao (1987). Power spectra of fluid velocities measured by laser Doppler velocimetry. Exper. in Fluids, 5, pp. 17-28. Benedict, L. H., H. Nobach and C. Tropea (2000). Estimation of turbulent velocity spectra from Laser Doppler data. Measurement Sci. Technol., 11, pp. 1089-1104. Bos, R., S. de Waele and P. M. T. Broersen (2002).
854
CONTINUOUS-TIME AUTOREGRESSIVE SPECTRAL ANALYSIS FOR IRREGULARLY SAMPLED DATA Piet M.T. Broersen Department of Multi Scale Physics, Delft University of Technology
Abstract: A continuous-time autoregressive spectral estimator is introduced that applies the principles of a discrete-time automatic equidistant missing data algorithm to unevenly spaced data. This time series estimator approximates the irregular data by a number of equidistantly resampled missing data sets, with a special nearest neighbor method. The ARMAsel-irreg algorithm estimates and automatically selects a discrete-time AR model from a number of candidates. This selected model often has a number of spurious high frequency poles, which are incompatible with the continuous character of the irregularly sampled signal. Those spurious poles can be eliminated, by transforming only the poles of the discrete time model with a positive real part to matching continuous-time poles. The estimated continuous-time spectra can be accurate at frequencies much higher than the mean data rate. Copyright © 2006 IFAC. Keywords: autoregressive model, nearest neighbor resampling, slotting, spectral estimation, time series analysis, uneven sampling, order selection, parametric model.
1. INTRODUCTION
normalization reduces the variance of the estimated autocovariance function. Fuzzy slotting produces a smoother autocovariance by distributing products over multiple time slots. Autocovariances estimated by slotting techniques are not guaranteed to be positive semi-definite. This results in spectra that can become negative at frequencies where the power is weak. Slotting methods require very large data sets.
Meteorological data, astronomical data or turbulence data obtained by Laser-Doppler anemometry are often irregularly sampled, due to the nature of the observation system. This has the advantage that the highest frequency that can be estimated is higher than half the mean data rate, which is the upper limit for equidistant observations. Many estimation techniques exist for unevenly spaced data. Lomb (1976) and Scargle (1989) described a Fourier estimator. Bos, et al. (2002) compare several variations of slotting or resampling methods. Larsson and Söderström (2002) replace the derivative operator in the continuous-time model by a discrete time approximation and apply it to a second order process. Lahalle, et al. (2004) estimate continuous time ARMA models, requiring the explicit use of a model for the irregular sampling instants. Larsson and Larsson (2002) give the Cramer-Rao lower bound for continuous-time ARMA models.
Resampling techniques reconstruct a signal at equal time intervals. After resampling, the equidistant data can be analyzed using the periodogram or time series models. Spectral estimates at higher frequencies will be severely biased. Adrian and Yao (1987) described Sample and Hold reconstruction as low-pass filtering followed by adding noise. These effects can in theory be eliminated using the refined Sample and Hold estimator (see Benedict, et al., 2000). In practice, all spectral details smaller than the bias are lost. Nearest Neighbor resampling has similar characteristics. The spectra are strongly biased for frequencies higher than about 20 % of the mean data rate. The noise and filtering effects of equidistant resampling set limits to the achievable accuracy of resampling methods. This precludes the accurate estimation of spectra at higher frequencies where the resampling noise blurs details,
Benedict, et al. (2000) give an extensive survey of methods for irregular data. Slotting methods estimate an equidistant autocovariance function from irregular samples. Slotting has been refined with normalization and fuzzy slotting. Local 849
smaller than the bias and hides spectral slopes.
any slot. Therefore, multi shift resampling is used, where M different equidistant missing data signals are extracted from one irregular data set and all small slots of width w are connected in time.
Bos, et al. (2002) introduced an estimator that can be perceived as searching for uninterrupted sequences of data that are almost equidistant. The contiguous sequences of different lengths can be analyzed with an irregular version of the Burg (1967) algorithm for segments. Slotted nearest neighbor resampling creates an equidistantly sampled signal by replacing an irregular observation time instant by the nearest equidistant resampling grid point, but only if it is within half the slot width. The slot width is equal to the resampling time Tr. Some grid points may remain empty. A smaller slot reduces the bias. The bias of slotted nearest neighbor is very much smaller than the bias without slotting. The reason is that one original irregular observation can never appear at multiple resampled time instants if slotting is applied. A disadvantage of this slotted Burg resampling method is that still very large data sets are required to obtain some uninterrupted equidistant sequences of sufficient length for the irregular Burg algorithm. It turned out that a non-linear maximum likelihood algorithm for missing data, developed by Jones (1980), could give a better solution, also if a limited amount of data is available.
This paper studies the possibility of automatic AR order selection for irregular data. It also describes how the transformation of the discrete-time model, with the automatically selected order, to a continuous-time model is favorable for the spectral accuracy. By simply eliminating poles of the selected estimated discrete-time model that cannot possibly be produced by poles of the continuous-time data, an accurate autoregressive spectrum is obtained. 2. TIME SERIES MODELS A discrete-time AR(p) model can be written as (Priestley, 1981) xn a1 xn 1 a p xn p
(1)
Hn ,
where Hn is a purely random process of independent identically distributed stochastic variables with zero mean and variance VH2. Assume that data represent a stationary stochastic process. The power spectral density h(Z ) of an AR(p) model is completely determined by the parameters in (1) together with the variance VH2 and is given by:
To study the possibilities of discrete-time models for irregular data, a separate study has first been made to the simpler problem of equidistant observations with missing data. In that case, Jones (1980) described an efficient method to calculate the true likelihood. A survey of existing missing data methods and a numerically stable version of the maximum likelihood algorithm for autoregressive models of missing data problems has been given by Broersen, et al. (2004a, 2004b). The performance of that automatic time series algorithm outperforms all other methods for missing data problems. The best method for missing data will be applied to irregular data.
h(Z )
V H2 1 2S A e jZ i p
2
V H2 2S
1 p
2
.
(2)
1 ¦ ai e jZ i i 1
Furthermore, the invariance property of maximum likelihood theory shows that the transformation (2) of efficiently estimated time series parameters also gives an efficient estimate of the power spectrum. Broersen (2002) developed a computer program to estimate the parameters and to select the model order with no other input than the data. It requires no user interaction to determine the best model for a given set of random observations.
Broersen and Bos (2006) showed that the continuous time maximum likelihood approach of Jones (1981) mostly fails to converge for irregular observations. If more than 2 parameters have to be estimated, the numerical surface of the continuous likelihood is too rough for global convergence, due to local maxima. Modifications are required to apply the automatic algorithm for equidistant missing data of Broersen, et al. (2004a,b) to irregularly sampled data. That needs multi shift slotted nearest neighbor resampling and uses a slot width that is only a fraction of the resampling time. This produces several equidistant data sets, with slightly different starting points that are shifted over the slot width. Time instants are fixed to an equidistant resampling grid where the original irregular sampling instant is never further away from the grid point than half a slot width. The resampling grid time Tr determines the discrete-time frequency range. By taking a slot width w that is only a fraction 1/M of the resampling time Tr, the bias is further reduced. This, however, gives disjoint intervals where many irregular times ti are not within
A similar automatic, approximate maximum likelihood, ARMAsel-mis program has been developed for the equidistant missing data problem. The accuracy of the ARMAsel-mis spectra is better than the spectra obtained with many other methods from the literature, in simulations with missing data. Broersen, et al. (2004b) have given examples where the estimation of time series models in missing data problems was efficient, meaning that the accuracy of the resulting model approached the limit of the achievable accuracy. The application of that discretetime missing data program to continuous-time irregular data will be studied here. That is necessary because the irregular maximum likelihood method of Jones (1980) mostly fails to converge and no other reliable automatic spectral estimation methods for finite or small samples of irregular data are available. The missing data AR order selection criterion GIC ( p )
850
LH D p,
(3)
uses the minimum negative log likelihood LH with a penalty factor D that depends on the missing fraction, 3 for less than 25 % missing, 5 for more than 75 % missing and 4 in the range between. T0 /MTr approximately gives the remaining fraction where T0 equals 1/f0, which in turn is the mean data rate.
the integral of the squared difference between the logarithms of spectra can be computed for an arbitrary frequency range 1/2Ta :
A continuous-time AR(p) model can be written as (Priestley, 1981) (4) s p D1s p1 D p x(t ) H (t ) , sx(t ) dx dt where s denotes the differentiation operator and H(t) is white noise with as correlation a G function and as spectrum the constant V 2. The power spectral density h(Z ) of an AR(p) model is given by the parameters :
The hat indicates a spectral estimate. By limiting the integration to the definition area 1/2Tr of the discrete spectrum, it is possible to attribute a number to the accuracy of discrete-time or continuous-time approximations.
h(Z )
V2 2
p
¦D
i
e
, D0
1, f Z f .
SD
NTa 2S
S / Ta
³ S
/ Ta
2
^
`
ªln ^h(Z )` ln hˆ(Z ) º dZ . ¬ ¼
(6)
3. ARMASEL FOR IRREGULAR DATA In the first application of the ARMAsel-irreg algorithm, Broersen and Bos (2006) reported good results for known low order models. Automatic order selection sometimes failed because of the peculiar behavior of the discrete-time likelihood function as a function of the model order. Often, spurious high frequency peaks seem statistically significant. That is illustrated in Fig. 1 for an example with AR(1) data.
(5)
jZ ( p i )
i 0
The discrete-time spectrum in (2) covers the frequency range Z between -S and S; the continuoustime spectrum has Z between -f and f. Transformations between models must always be a compromise of properties. Franklin, et al. (1998) describe transformations from a control point of view. The bilinear or Tustin transformation is often a good choice for the conservation of filter properties. A disadvantage is that the frequency scale is warped. The zero order hold equivalents can give the same outputs for continuous-time and discrete time systems at the sampling instants. Söderström (1991) gives several solutions for the computation of stochastic continuous-time models from discretetime ARMA models. He showed that negative real discrete-time poles and also some more specified areas of zero positions do not give any equivalent continuous-time solution. Those stationary and invertible ARMA models cannot possibly be obtained by sampling a continuous-time process.
Fig.1. Spectrum of discrete-time AR models, shifted downwards with a factor 10 for each new order. 8 AR models are estimated from N=100 irregular samples of the true continuous AR(1) process with the pole at - 0.1S, resampled with Tr = 1/2f0, where f0 =1 is the mean data rate; w = Tr. In this example, the AR(3) model was selected.
For spectral analysis, the conservation of the location of spectral peaks and valleys is the most important. That can be obtained with a technique called pole and zero matching. A pole at s = D in (4) will give a pole at exp(DTr) in (1). A pole at s = D + iE will give a pole at Uexp(iM) in (1), with U = exp(DTr) and M = ETr, (Franklin, et al., 1998). Stationary or stable processes must have only poles with a negative real part in the continuous s-domain or a radius less than 1 in discrete time models. The same matching rules can be used for the inverse transformation of discrete-time to continuous-time models. All discrete poles within the unit circle have a stable transform, except negative real poles for which no sensible matched pole transform can be defined. AR processes have no zeros, otherwise the same rules would be used for the matched transform of zeros.
The order 3 was selected in Fig. 1, although the best approximation was given by the AR(1) model. It is curious that the AR(3) model has a strong peak at the end of the frequency interval. This belongs to a negative real pole, which cannot belong to a sampled continuous-time process. All the higher order models show spurious details. Strange peaks can mostly be found in the range between 1/4 Tr to 1/2Tr. What is worse, models with such ‘impossible’ peaks at high frequencies are often selected by GIC(p) of (3).
Accuracy measures have been defined by de Waele and Broersen (2000) that can compare the quality of the discrete-time model with true or with aliased continuous spectra. With the spectral distortion SD,
Simulations with higher order generating processes have a similar behavior. Estimated low order models may have acceptable spectra. But at and above the true order, and often already below, estimated
851
models often have spurious peaks that seem to be statistically significant. Most peaks appear at the highest discrete frequency 1/2Tr, many in the last half of the frequency range, above 1/4Tr, and only very few are found for lower frequencies. Order selection criteria often select one of the discrete-time models with a peak at 1/2Tr, a model without a continuous time equivalent. It has been verified that this agrees with the minimum of the log likelihood in those examples. In other words, the likelihood decreases with a statistically very significant amount by adding a negative real pole, giving a spectral peak at the frequency 1/2Tr. No value of the penalty factor D in (3) would give a good performance of the order selection criterion. An analysis of the bias of ARMAsel-irreg by Broersen and Bos (2006) showed that the frequency range between 1/4 Tr to 1/2Tr suffers from serious bias in examples with strong spectral slopes at those frequencies. That frequency range is very sensitive to bias and to spurious peaks.
4. PERFORMANCE OF THE NEW ESTIMATOR Simulations with a known spectrum are a first step in testing new algorithms. Test data was generated with N irregular continuous-time data points of a known AR process. The resulting irregular data has time intervals that are Poisson distributed. In simulations, the true properties of the data are known. Hence, the quality of estimated results can be established with the SD of (6). The vertical scales in the spectral plots are sometimes shifted, to increase the visibility.
An example can be given where a missing data problem with a pole at a1 could have a minimum of the likelihood for a pole at – a1. That can happen for an AR(1) process with parameter a1 if the data are missing with a certain pattern. Suppose that all odd observations are missing, only even time instants remain, for T = 1. The autocorrelation function of this process is the series (-a1)k and that is identical for the parameters + a1 or – a1, for all even values of k. The convergence of the optimization algorithm for the AR(1) parameter depends on the starting values then, because the likelihood is symmetric. That has been verified. The same regular pattern of only even samples could theoretically also be found for irregular sampling. It gives at least an explanation that it is possible to find those strange, ‘impossible’, peaks in an algorithm without programming errors.
Fig. 2. True spectrum, discrete-time ARMAsel-irreg and two continuous-time estimates for a true AR(2) spectrum. Tr = 1/4f0, with f0 = 1 and slot width w = Tr /2. N = 250. ARMAsel-irreg selected AR(4). Fig. 2 gives spectra for 250 irregular observations of an AR(2) process with 1 peak. The ARMAsel-irreg program selected the discrete-time AR(4) model from candidates AR models until order 6, with a log likelihood value LH equal to - 152.8. The estimated discrete-time AR(2) model had LH= - 122.3. That spectrum fitted much better than the selected AR(4) model. However, after matched transformation to continuous-time, the ARMAsel continuous Matched AR(4) model had a better spectral estimate than the transformed AR(2) model. The matched continuous model eliminates all discrete poles with negative real parts, which agrees with frequencies greater than Tr/4. The bilinear continuous transformed spectrum keeps the wrong peak of the ARMAsel-irreg AR(4) model and shifts it to a new frequency in Fig. 2. This shifting property as well as the conservation of the peak with a frequency above Tr /4 makes the bilinear transformation less adequate for this application. Repeating this simulation 10 times always gave a selected AR(3) or AR(4) model, with about the same accuracy for the final ARMAsel continuous matched model. The true continuous process has a spectrum with an infinite frequency range, the estimated discrete-time model has a frequency range until 2f0 in Fig. 2. The range of the transformed continuous-time models is also infinite. However, as those models are derived from the discrete-time model, it is an open question how reliable high frequencies are. An
Two completely different algorithms for the log likelihood for missing data have been programmed independently by Broersen, et al. (2004a, 2004b): the exact likelihood as formulated by Jones (1980) and an approximation. Both algorithms lead to the same estimates for the likelihood and to the same spurious high frequency spectral peaks. Broersen and Bos (2006) also showed that the continuous likelihood has an undesirable character, a conclusion that was also supported by other authors. It still is an open question whether the regular even sampling pattern of the previous paragraph is a possible explanation for peculiar behavior in Fig. 1. No other explanation is available yet. A practical solution for the order selection problem for irregular observations will be investigated by switching to continuous time models. First, ARMAsel-irreg selects the discrete-time model with the minimum of GIC(p) of (3). Afterwards, all poles of the selected model with frequencies above 1/4Tr will be eliminated before a matched transform is applied. By choosing the value of Tr small enough, it can be made sure that no true spectral peaks are in that frequency range. In this way, the frequency range that is sensitive to bias and to spurious spectral peaks is eliminated in the continuous-time model.
852
Fig. 4 gives the result for an AR(4) process with two strong spectral peaks. All discrete-time higher order models showed a strong peak at the end of their frequency range. In this example the AR(4) model would also give a good spectrum. In another run with 10000 observations, the estimated AR(4) model found only the second true peak, and an extra peak at the end of the frequency range. It missed the first true peak. Obviously, the spurious high frequency peak was more significant in the likelihood. In that run, the selected model order was much better than the true order. That often has been seen in individual simulation runs. Knowing the true order mostly has no advantage for irregular data in our simulations. Using selected orders is generally more accurate.
intermediate range is used for the visual representation, to indicate the differences. The perfect match of the estimated Matched spectrum at the high frequencies is caused by the agreement of the number of true and estimated poles. Each pole gives a high frequency spectral asymptote of –20 dB per frequency decade in continuous modeling.
It has been tried to use the parameters of the selected continuous matched ARMAsel model as starting values for a further search of a continuous maximum likelihood model. The very irregular surface of the continuous likelihood function, as reported by Broersen and Bos (2006) prevented this approach. Sometimes, the estimated model did not improve but remained the same. It could also happen that the continuous likelihood converged to a much poorer model than the starting position. It never occurred in the simulations that the AR(4) example converged to a much better model. This continuous maximum likelihood post processing could only be used fruitfully for AR(1) and AR(2) examples, where also the approach of this paper did not give any problem.
Fig. 3. True spectrum, ARMAsel-irreg and two identical continuous matched estimates for AR(3) data. Tr = 1/2f0 with f0 = 1 and slot width w = Tr /2 , N = 1000. ARMAsel-irreg selected AR(3). Fig. 3 gives the result for a smooth AR(3) process with the three real poles at -0.1S, -0.4S and -S. AR(3) was selected with ARMAsel-irreg and higher order models showed only a moderate tendency for spurious details. In general, examples with strong details in the spectra have a tendency to overestimate the order and to have more spurious peaks for the same sample size. In this example, the selected and the true order model coincided because AR(3) was selected. Remind the vertical shift of other spectra.
Fig. 5. True, discrete-time ARMAsel-irreg spectrum, and a continuous matched estimate for a true AR(5) spectrum. Tr = 1/2f0 with f0 = 1 and slot width w = Tr /2. N=2000. ARMAsel-irreg selected AR(3). Fig. 5 has a rather smooth true spectrum with 5 real poles at -0.1S, -0.2S, -0.4S, -0.8S, and -1.6S. AR(3) was selected by ARMAsel-irreg. The estimated AR(5) model had one negative real pole, which was eliminated. The different orders can be recognized in the different asymptotic slopes, -60, -80 and –100 dB for AR(3), AR(4) and AR(5), respectively. The performance of this estimator and the spurious peaks for greater sample sizes still have to be investigated.
Fig. 4. True, discrete-time ARMAsel-irreg spectrum, and a continuous matched estimate for a true AR(4) process. Tr = 1/4f0 with f0 = 1 and slot width w = Tr/2. N = 1000. ARMAsel-irreg selected AR(6).
853
Autoregressive spectral estimation by application of the Burg algorithm to irregularly sampled data. IEEE Trans. on Instrumentation and Measurement, 51, pp. 1289-1294. Broersen, P.M.T. (2002). Automatic spectral analysis with time series models. IEEE Trans. on Instrument. and Measurement, 51, pp. 211-216. ARMASA Matlab toolbox, freely available at www.mathworks.com/matlabcentral/fileexchange.
Broersen, P.M.T. and R. Bos (2006). Estimation of time series models from irregularly spaced data. IEEE Trans. on Instrumentation and Measurement, offered for publication. Broersen, P.M.T., S. de Waele and R. Bos (2004a). Application of autoregressive spectral analysis to missing data problems. IEEE Trans. on Instrument. and Measurement, 53, pp. 981-986. Broersen, P.M.T., S. de Waele and R. Bos (2004b). Autoregressive spectral analysis when observations are missing. Automatica, 40, pp. 1495-1504. Burg, J. P. (1967). Maximum entropy spectral analysis. Proc 37 th Meeting Soc. Of Exploration Geophysicists, 6 pp. Franklin, G.F., J.D. Powell and M. Workman (1998). Digital Control of Dynamic Systems. Addison Wesley, Menlo Park, CA, USA. Jones, R.H. (1980). Maximum likelihood fitting of ARMA models to time series with missing observations. Technometrics, 22, pp. 389-395. Jones, R.H. (1981). Fitting a continuous time autoregression to discrete data. Applied Time Series Analysis II, (Ed. D.F. Findley), pp. 651682. Lahalle, E., G. Fleury and A. Rivoira (2004). Continuous ARMA spectral estimation from irregularly sampled observations. Proc. IEEE/IMTC Conf., Como, pp. 923-927. Larsson E.K. and T. Söderström (2002). Identification of continuous-time AR processes from unevenly sampled data. Automatica, 38, pp. 709-718. Larsson, E.K. and E.G. Larsson (2002). The CRB for parameter estimation in irregularly sampled continuous-time ARMA systems. IEEE Signal Processing Letters, 11, pp. 197-200. Lomb, N.R. (1976). Least squares frequency analysis of unequally spaced data. Astrophysics and Space Science, 39, pp. 447-462. Priestley, M. B. (1981). Spectral Analysis and Time Series. Academic Press, London, U.K. Scargle J.D. (1989). Studies in astronomical time series analysis III. Fourier transforms, autocorrelation functions, and cross-correlation functions of unevenly spaced data. The Astrophysics Journal, 343, pp. 874-887. Söderström, T. (1991). Computing stochastic continuous-time models from ARMA models. Int. J. Control, 53, pp.1311-1326. Waele S. de and P. M. T. Broersen (2000). Error measures for resampled irregular data. IEEE Trans. on Instrumentation and Measurement, 49, pp. 216-222.
Fig. 6. True, discrete-time ARMAsel-irreg spectrum, and two continuous estimates for a true AR(6) spectrum. Tr = 1/4f0 with f0 = 1 and slot width w = Tr /2 , N = 2000. ARMAsel-irreg selected AR(10). The example in Fig. 6 demonstrates that sometimes the spurious peaks are estimated in AR(6) before the true peaks. The elimination of the discrete-time poles with negative real parts from AR(10) is very favorable here. AR(10) was the highest candidate.
5. CONCLUSIONS A new estimator is introduced that first fits discretetime AR models to multi shift slotted nearest neighbor resampled segments from irregular data. The ARMAsel-irreg algorithm automatically selects the AR order. Afterwards, the reliable poles of the selected discrete-time AR model are transformed to matched continuous-time poles. Only estimated poles with a positive real part are used. All poles with negative real parts are eliminated, because they mostly belong to spurious spectral details. In simulations with few irregular data, the results are much better than what can be obtained with any other known existing spectral estimation technique. Multi shift slotted nearest neighbor resampling will give accurate spectra, with a small bias, if the dynamic spectral range is limited. For a large dynamic range, a small slot width reduces the expected bias. The elimination of continuous-time poles with a frequency above 1/4Tr further reduces the bias. In practice, Tr can be chosen small enough to make sure that no true process poles will be eliminated. REFERENCES Adrian R. J. and C. S. Yao (1987). Power spectra of fluid velocities measured by laser Doppler velocimetry. Exper. in Fluids, 5, pp. 17-28. Benedict, L. H., H. Nobach and C. Tropea (2000). Estimation of turbulent velocity spectra from Laser Doppler data. Measurement Sci. Technol., 11, pp. 1089-1104. Bos, R., S. de Waele and P. M. T. Broersen (2002).
854