- Email: [email protected]
Correlation estimation and time-series modeling for nonstationary processes
Signal Processing 15 (1988) 31-41 North-Holland
31
C O R R E L A T I O N E S T I M A T I O N A N D TIME-SERIES M O D E L I N G FOR N O N S T A T I O N A R Y P R O C E S S E S * William A. G A R D N E R (Member EURASIP) Signal and Image Processing Laboratory, Department of Electrical Engineering and Computer Science, University of California, Davis, CA 95616, U.S.A. Received 18 February 1987 Revised 9 November 1987
Abstract. Various types of statistical signal processing require the estimation of time-variant correlation functions and the fitting of time-variant models for nonstationary processes. In this paper, the situations for which nonstationary probabilistic correlation functions can be accurately estimated from a single sample path of a stochastic process are delineated. These include only (i) known form of nonstationarity, (ii) periodic or almost periodic nonstationarity, and (iii) slowly fluctuating nonstationarity. Known methods of estimation for each of these situations are reviewed within the unifying framework of orthogonal series expansions. The fact that estimators based on orthogonal series expansions for nonstationary correlation functions other than (i)-(iii) cannot be guaranteed to be accurate is established. Ramifications for time-variant autoregressive model fitting are discussed. Zasammenfassung. Verschiedene Aufgabenbereiche der Signalverarbeitung erfordern die Sch~itzung zeitvarianter Korrelationsfunktionen und die zeitvariante Modellierung nichtstation~irer Prozesse. lm folgenden werden die Situationen dargelegt, fiir die nicht-station~ire Korrelationsfunktionen aus einem einzigen Repr~isentanten eines Prozesses gesch~itzt werden k/~nnen. ErfaJ3t werden kSnnen nur F~ille mit (i) bekannter Form der Instationarit§t, (ii) periodischer oder fast-periodischer Instationaritat und (iii) langsam ver/inderlicher Instationarit~it. Die bekannten Sch~itzverfahren fiir jede dieser Gegebenheiten werden in einem einheitlichen Rahmen, dem der Orthogonalreihen-Entwicklung, zusammengestellt. Die Grenzen der Nutzung von Reihenentwicklungen bei der Sch~itzung nicht-station~irer Korrelationsfunktionen werden aufgezeigt. VariationsmSglichkeiten fiir die zeitvariante AR-Modellierung werden diskutiert. Rrsumr. De nombreuses operations de traitement du signal requi~rent l'estimation de fonctions de corrrlation dependant du temps et la construction de modules variant avec le temps pour les processus nonstationnaires. Dans cet article, les cas pour lesquels les fonctions probabilistes de corrrlation nonstationnaires peuvent ~tre estimres de mani~re pr6cise a partir d'une seule trajectoire d'un processus stochastique sont repertorirs. Ces cas incluent seulement: (i) les formes connues de stationnaritr, (ii) les processus p6riodiques ou presque prriodiques, et (iii) les processus nonstationnaires ~ evolution lente. Les m&hodes connues d'estimation pour chacune de ces situations sont reexaminres du point de vue grnrral des reprrsentations en srries orthogonales. Le fait qu'il n'est pas possible de guarantir la prrcision des estimres sous forme de srries orthogonales pour les fonctions de corrrlation nonstationnaires autres que du type (i)-(iii) est 6tabli. Les cons6quences de ce rrsultat pour la construction de modules autorrgressifs drpendant du temps sont examinres. Keywords. Nonstationary correlation, correlation estimation, time-series modeling, nonstationary time-series.
1. Introduction ~ Various types of signal processing require the estimation of time-variant correlation functions and the fitting of time-variant models for nonstationary processes. The most straightforward * This work was supported in part by a grant from ESL, Inc., Sunnyvale. CA.
way to estimate the correlation function of a nonstationary stochastic process using one sample path x(t) over a time interval [0, T] is to use a sliding window average of the lag product l The research reported by W. Music in "Recursive prediction and estimation of periodic nonstationary signals", 1987 IEEE Region 5 Conference, pp. 139-142, duplicates that presented herein, but was mistakenly not attributed to the author of the present paper.
0165-1684/88/$3.50 © 1988, Elsevier Science Publishers B.V. (North-Holland)
32
W..A. Gardner / Nonstationary correlation estimation
x ( t ) x ( t - r ) . When possible, the window must be narrow enough to resolve the nonstationarity but wide enough to provide adequate averaging. An alternative approach proposed by Snyder [20] uses an orthogonal series expansion in t of the lagproduct time-series x ( t ) x ( t - r ) over the entire interval [0, T], with r treated as a parameter. Boyles and Gardner [1] showed that the orthogonal series expansion method yields consistent estimators (mean squared error approaches zero as T-~ oo) for periodic and almost periodic nonstationarity (cyclostationarity) when the complex sinusoids with appropriate frequencies are used as the basis functions, provided that the process exhibits appropriate ergodic properties (cycloergodicity). 2 But the utility of this approach for more general nonstationarity has not been established. This orthogonal series expansion approach can be applied to autoregressive modeling of nonstationary processes by simply substituting the nonstationary correlation estimates for the true nonstationary correlations into the probabilistic normal equations whose solution provides a minimum mean squared error autoregressive model. A related alternative approach uses orthogonal series expansions of the time-variant autoregressive model parameters and then derives empirical normal equations whose solutions provide the least squares estimates of the expansion coefficients. This approach has been proposed and studied by many authors including Jones and Brelsford [9], Subba Rao [21], Liporace [ 11], Hall et al. [7], Parzen and Pagano [16], Hasselman and Barnett [8], Grenier [6], and Sharman and Friedlander [18]. It can be argued that the attainable performance of this alternative approach is tied to that of the former approach. Thus, insight into the general utility of orthogonal series expansions for both nonstationary correlation estimation and nonstationary autoregressive modeling can be gained by studying the former problem which is the more analytically tractable of the two. 2An almostperiodicfunction is a sum of periodic functions with incommensurateperiods. Signal Processing
With the wide-spread applicability of any method for obtaining good estimates of nonstationary correlation functions as motivation, the purposes of this paper are (i) to delineate the types of nonstationarity for which consistent or at least accurate estimates of nonstationary correlation functions can be obtained and to explain how to obtain them, and (ii) to establish that orthogonal series expansions cannot provide estimates of nonstationary correlation functions that are guaranteed to be accurate for any other types of nonstationarity. The remainder of this introduction provides a summary of the current status of the problem of estimation of nonstationary correlation. The conventional method for estimating a correlation function, say Rxy(7. ) a E{x(t)y(t - ~')},
(1)
for a pair of jointly wide-sense stationary stochastic processes x(t) and y(t) using one sample path from each of x(t) and y(t) is to time-average the lag product:
a if ofx ( t ) y ( t - ~') dt.
Rxy(7) =-~
(2)
If x(t) and y(t) possess appropriate ergodic properties (see [4, pp. 150-153] and Section 2), then we have the desired limit A
lim Rxy(r) = Rxy(r)
T~oo
(3)
in some probabilistic sense such as in mean square or with probability one (depending on the particular ergodic property). It has been shown by Parzen [15] that (3) also holds if x(t) and y(t) are only asymptotically jointly wide-sense stationary. In fact, even if x(t) and y(t) are generally nonstationary, but are asymptotically mean stationary in the sense that the time-average Rxy(r) & lim 1 S~oo S
Rxy(t, r) dt
(4)
of the nonstationary correlation function
Rx~,(t, r) a E { x ( t ) y ( t - r)}
(5)
33
W.A. Gardner / Nonstationary correlation estimation
exists and is not identically zero, then (3) holds provided that x(t) and y(t) possess appropriate ergodic properties (see [10], [4, pp. 150-159], and Section 2). More recently, it has been shown by Boyles and Gardner [ 1] that by using an asymptotically orthonormal series expansion the nonstationary correlation function (5) itself (not just its timeaverage) can be consistently estimated (mean square estimation error approaches zero for every t, or almost every t, as averaging time T approaches infinity) provided that the nonstationarity is either cyclostationary or almost cyclostationary (Rxy(t, z) is either periodic or almost periodic--in the mathematical sense---in t) and provided that x(t) and y(t) possess appropriate ergodic properties called cycloergodicity (see also [4, pp. 345-349] and Section 2). It has been claimed by Marmarelis [13, 14] that the orthogonal series method proposed by Synder [20] can provide consistent estimators for more general nonstationary correlation functions, However it is shown here in Section 3 that this is so only if the form of nonstationarity is known in advance.
2. The orthonormal expansion method In order to analyze the consistency of the estimator to be described it must be assumed that the stochastic process models are sufficiently well behaved that the estimator actually exists. This can be accomplished by requiring that the processes x(t) and y(t) be measurable so that almost all sample functions of the lag-product x ( t ) y ( t - T ) are appropriately integrable. Also, in order to prove that the mean squared error of an estimator obtained from an orthonormal expansion approaches zero for every value of t as the integration time T approaches infinity, it must be assumed that R~y(t, r) is a continuous function of t; otherwise only the L 2 norm of the error can be expected to approach zero or the error can be expected to approach zero for every value of t except for a set of zero measure. Of course, for any physical phenomenon for which the correlation estimation
problem makes sense, there will exist an appropriate stochastic process model. Consequently, these technical aspects of the problem are minimized in the treatment presented here since the primary objective is to present arguments and results that are directly useful to the practitioner. Motivated by the classical method of orthonormal expansions of a function, the following method for estimation of a nonstationary correlation function is considered. Let {bm(t)}-MM be a set of orthonormal functions on the interval [0, T], llo r
bm(t)b*(t)dt=
{1,
m=n,
O, m # n.
(6)
The proposed estimator for Rxy(t, r) is given in terms of {b,,(t)}_MM by M
J~xy(/, 7") & ~, t~m(z)bm(t),
(7a)
m=-~ where a,,(T) =-~ ^ ~ llor x(t)y(t-z)b*(t)dt.
(7b)
In the special case for which M = 0 and bo(t) = 1, this method reduces to the conventional method for estimating a stationary correlation function Rxy(t, r)=- Rxy(z), or the time-averaged value,
ao(z)&~l for Rxy(t, ~') dt,
(8)
of a nonstationary correlation function. It is wellknown that the method produces consistent estimators (mean squared error approaches zero as T - * ~ ) in these two cases provided that the stochastic processes x(t) and y(t) exhibit the appropriate ergodic property (see [2] and [4, pp. 150-159]), namely (14) for b,,(t) =- 1. In the special case for which {bm(t)}_MM are sine waves, bm(t)=e i2~m',
ot-,,,=-Otm,
(9)
this method reduces to that studied by Boyles and Gardner [1], where it is shown that it produces consistent estimators provided that: Vol. 15, No. l, July 1988
34
W . A . G a r d n e r / N o n s t a t i o n a r y correlation e s t i m a t i o n
(i) the correlation function R x y ( t , 7") is either periodic ({am}_MM are all integer multiples of a single fundamental frequency) or almost periodic with only M harmonics, M
Rxy(t, 7,)= Y~ am(7.)e i2~ nl,, rn=
M
it can be shown by methods analogous to those used by Boyles and Gardner [1] that the estimator (7) is consistent provided that x(t) and y(t) possess appropriate ergodic properties. This is accomplished by showing that the bias is asymptotically zero, lim E{/~xy(t, ¢¢)}-Rxy(t, 7-)---0,
(12)
a_m(7,)= a*(7.), (lOa) for almost all t and 7-, and that the variance of each one of the coefficient estimators in (7) is asymptotically zero,
where am(~') ~ lim 1 S~oo S
Rxy(t, 7,) e
~ at,
i2 . . . .
(10b)
lim var{~m(7,)}-= 0,
T~oc5
(ii) x(t) and y(t) possess appropriate cycloergodic properties, namely (14) with bin(t) given by (9). In this case, the set {bm(t)}~_Mis not orthonormal unless {a,,}-M M are all integer multiples of a single fundamental frequency and T is an integer multiple of the reciprocal of this fundamental frequency. Nevertheless, the set {bm(t)}_~M is asymptotically (T-->oo) orthonormal and the estimator (7) is indeed consistent as shown by Boyles and Gardner [1] (see also [4, pp. 345-349]). In fact, it is shown by Boyles and Gardner [1] that even if M in (10) is not finite, the limit of the estimators (7) as T ~ o o and M ~ o o is consistent if (14) is satisfied. (For related work on consistent estimation of the Fourier transform of R x y ( t , • ), which is the instantaneous spectrum, see [5].) In the special case for which the correlation function is given by the finite separable form M
Rxy(t, 7.)= ~
a,,(7,)bm(t), O<~t<~,
m=-M
provided that property
x(t) and y(t) possess the ergodic
lfo fo
lim--~
Rz~,(t, u)b*(t)bm(u) dt du
-- 0, where
(13)
(14a)
z~(t) is the centered lag-product process
z~(t) a=x(t)y(t-7,)-E{x(t)y(t-7-)}. (14b) For example, x(t) and y(t) can be jointly Gaussian with correlation functions that asymptotically decay to zero as the lag 7 approaches infinity. (This is demonstrated in the next section.) Of course, the fact that {bm(t)}M_M in (7) and (11) must be one and the same implies that the form of nonstationarity must be known before the estimate (7) can be formed. For example, if
x(t) = bo(t)y(t) and y(t) is stationary, then (1 la) holds with M --- 0, and if bo(t) is persistent in the sense that
(lla) rlim ~ IT
where
for ]bo(t)[ 2 dt = 1,
then (llc) holds, and ( l l b ) becomes
Rxy(t, 7.)b*(t) dt, am(7")~-limTIo r r ~
(llb)
and the functions {bm(t)}_~M are asymptotically orthonormal lim
T~c~
T
lfo
Signal Processing
bm(t)b*(t)dt=
{1,
re=n,
0,
m ~a~ n ,
(llc)
ao(r) = Rry(Z). Furthermore, if y(t) exhibits mean square ergodicity of its autocorrelation (see [4, pp. 150-153]; e.g., y(t) could be Gaussian with no additive random sine wave components), then the estimator (7) with M = 0 is consistent.
W.A. Gardner / Nonstationary correlation estimation
Another example which better illustrates the constraint (11c) is
x(t) = bo(t)Zo(t) + bl(t)zl(t) + b_l(t)z_l(t), where Zo(t), zl(t), z_l(t), and y(t) are zero-mean jointly stationary ergodic Gaussian processes, bo(t)-=l, and bl(t) and b_l(t) are two +1 waveforms whose zero-crossings are specified by two independent samples of a Poisson point process. Then ( l l c ) as well as ( l l a ) and ( l l b ) are satisfied. This illustrates that a wide variety of appropriate {b,,(t)}ffM other than sine waves do indeed exist. The real limitation to the practical applicability of the required form (11) is that the basis functions {b,,(t)}ffM must be known in advance. More generally, if (11a) is only an approximation, but this approximation becomes exact in the limit as M -> oo, then the estimator (7) is consistent in the limit M--> oo, if (14) is satisfied. Or, in other words, for any arbitrarily small positive number e, there exists a positive integer M and a positive number T both sufficiently large to guarantee that the mean squared error of the estimator (7) is less than e for all t e [0, T] and all • such that t - ~"e [0, T]. This is proved for the special case of sinusoidal basis functions (9) by Boyles and Gardner [1]. This suggests that the estimator (7) might be made consistent for any unknown but persistent nonstationarity for which (14) is satisfied by letting M-> oo and using any complete orthonormal set on L2[0, T] (the space of square integrable functions on the interval [0, T]) that is parameterized by T so that (11c) is satisfied. The objectives of the next section are to show that, for general unknown nonstationarity, such an estimator is inconsistent and for finite T and M cannot be guaranteed to be accurate. However, before proceeding, it is explained that for nonstationarity that is neither of known form nor periodic (nor almost periodic), the only situation in which accurate estimates can be guaranteed is that where the nonstationarity is known to be sufficiently slow, typically called local stationarity or quasi-stationarity. It is shown in the next
35
section that in this case (7) is appropriate with only M = O and bo(t)=-1 (although M > 0 can be used as explained by Snyder [20]). Generally speaking, a pair of processes is said to be jointly locally stationary (see [19] and [20]) if the bandwidth B of their nonstationarity and the maximum width ro of their auto- and cross-correlation functions are related by Bro,~ 1.
(15)
The parameter To is a measure of the width of, say, Rxy(t, r) as a function of ~', and B is a measure of the spectral bandwidth of, say, Rxy(t, ~') as a function of t. If the averaging time T in the conventional estimate (2) is chosen to satisfy
T
T~'~o,
(16a, b)
then the mean squared error between /~xy(r) and Rxy(t, r) will be small for 0<~ t ~< T, provided that x(t) and y ( t - s ) are nearly statistically independent for Isl > In order to track the slow nonstationarity, the fixed window [0, T] can be made into a sliding window [t', t ' + T] to obtain the estimate A , allot Rxy(t,'r)=-~
x(t'+t)y(t'+t-z)dt. (17)
The condition (16a) insures a small bias due to time smearing effects, and condition (16b) insures a small variance. (See [20] and, for related work on estimation of the Fourier transform of Rxy(t, • ) for locally stationary processes, see [12, 2, 17 and 5, and references therein].)
3.
Limitations
The fundamental reason why the estimator (7) cannot be made consistent with M ~ oo (and, of course, T - ~ ) for any correlation function Rxy(., .r) that is contained in L 2 [ 0 , T ] for all T < oo is that there is no separable space containing all such functions s. That is, it is not true that there 3 However, there is such a space for all almost periodic functions whose sine wave frequencies form a countable set. Vol. 15, No. 1, July 1988
W.A. Gardner / Nonstationary correlation estimation
36
exists a countable basis {bm(t)}_~oo depending on the parameter T such that for any such function Rxy( ", 7.) and any arbitrarily small positive number e there exists a positive integer M sutticiently large that
Rxy(t, 7.)-
Y~ am(7.)bm(t)
(18)
m=-M
for almost every t c [0, T] for arbitrarily large T, where am(7.) is given by ( l l b ) . What is true is that for a given accuracy e over the interval [0, T], the number M of terms used must grow with T if the basis {bm(t)}_MM is chosen arbitrarily. Consequently, as T is increased in order to reduce the variance of each term (7b) in the estimator (7a), the bias will grow unless M is also increased. But as M is increased, the variance of the sum in (7a) increases. In fact, for any fixed T, as M ~ ce, we obtain from Bessel's equality for a complete orthonormal set {bm(t)}_~oothe limit
lim Rxy(t, 7.)=x(t)y(t--7.),
M~oo
(19)
for almost every t ~ [0, T], which has an unacceptably large variance although the bias is zero. At the opposite extreme, it is shown in this section that for any fixed M, as T ~ , the chances of obtaining a small bias are M out c~, which is zero, although the variance is zero in this case. One approach to quantifying this problem of trading off bias and variance by selection of M and T is described in this section, and is used to illustrate the fact that (7) cannot be guaranteed to provide accurate estimates for general unknown nonstationarity. In order to prove that M must grow with T, we restrict consideration to the subset of all nonstationary correlations for which the nonstationarity is nearly bandlimited to arbitrarily large but finite bandwidth B. We then establish that M must grow linearly with T. (Since this restriction does not impose any limitation to the practical applicability of the result, we do not pursue a more general and necessarily more abstract proof that M must grow at a certain rate with T.) Signal Processing
Although the form of nonstationarity may be unknown, the physics of any practical situation would enable one to set an upper bound on the highest spectral bandwidth B of practical interest. This can be done mathematically by selecting B such that, say, 99% of the energy of every practically feasible function Rxy(., 7.) on the time interval [0, T] in a given practical situation is contained within the spectral band [ - B , B]. The smallest generally useful subset of L2[0, T] for purposes of approximation is defined to be the range space of the operator on L2(-oo, ~ ) that bandlimits to [ - B , B] and then duration limits to [0, T]. Thus, for every Rxy(" , 7") contained in this smallest generally useful subset, there exists a function f~(.) in L2(-oo, co) such that Rxy( ", 7.) is the image of f~(. ) under this operator:
Rxy(t, 7.)=f~_ook(t,u)f~(u)du,
t~[0, r], (20a)
where k(t, given by
u) is the kernel of the operator and is
k(t,u)=g(t)h(t-u),
(20b)
where
O<~t<~ T,
(20c)
otherwise is the duration-limiting part and
h(t - u) -
sin[27.rB(t - u)]
,.(t-u)
(20d)
is the bandlimiting part. Since this operator is compact, its range space is a compact subset of L 2 and is therefore approximately finite-dimensional. That is, there exists a finite dimensional subspace of L2[0, T] such that every function Rxy( ", 7.) contained in this compact subset has, say, 99% of its energy contained in this subspace. Specifically, the orthogonal projection, denoted by Rxy(', 7.), of Rxy(', 7.) onto this subspace satisfies II/~xy(.,
r)-R~y(., 7.)11
IIRxy(', 7.)11
<0.01,
W.A. Gardner/ Nonstationarycorrelationestimation where
Rxv(., z),, a=[T I ~ lRxy(t, z)]2 dt]~
(21)
is the L2[0, T] norm. If this finite dimensional subspace is spanned by the orthonormal basis {b,, (t)} ffM (cf. (6)), then the orthogonal projection is given by M
Rxy(t, z)=
Y~ am(z)bm(t),
(22a)
rn=-M
where
l forRxy(t, z)b*(t) dt.
a,,('r) = ~
(22b)
With the preceding definition in hand, we have the following theorem. Theorem. The smallest generally useful subset of L2[0, T] for approximation of functions with band-
width parameter B has approximate dimensionality D given by D = 4BT. The p r o o f of this theorem follows directly from the classical results (from the theory of functional analysis) on optimum finite-dimensional approximation of compact operators. A practically oriented account is given by Franks [3]. As a consequence of this theorem, the approximation M
E am(T)bm(t) ~Rxy(t, r) m ~ - M
oo
= E
am(z)b,~(t),
rtl = - c x 3
O<-t<~ T, (23) using an arbitrary basis {b,.}_~ for L2[0, T], cannot be expected to be accurate unless 2 M + 1 1>4BT (say M >12BT). In fact, even with M >!2BT (23) will not be an accurate approximation unless an appropriate basis (e.g. each b,.(t) approximately bandlimited to a spectral bandwidth of at most B) is used. To be most favorable to the proposed estimator (7), let us assume that the basis {bin} is
37
appropriate so that M need not be any larger than 2BT. If M is less than 2BT then the fact that a close approximation in (23) could be only a lucky chance insures that an estimate of the form (7) cannot be a close approximation to Rxy(t, z) except as a lucky chance. In fact one could say that the chances of (7) being capable of providing a good estimate are M out of 2BT. This chancy estimation is not of practical value since there is in general no way to know whether you have been lucky or not. Moreover, the chances of being lucky with a fixed M (which is all that can be used in practice) become smaller and smaller as T is increased, as explained in the following. However, before proceeding it is explained that if it is known that R~y(., z) is periodic with period To, then by choosing b,,(t) = exp(i2~rmt/To), (23) will be a close approximation for M=2BTo, regardless of how large T is. Thus, by choosing M = 2BTo to obtain a small bias, and choosing T sufficiently large to render the variance of the estimator (7a) small (assuming the ergodic property (14)), the estimator can be guaranteed to provide accurate estimates. A similar argument applies when Rxy(', r) is almost periodic. Also, when x(t) and y(t) are jointly locally stationary, then we can choose T to satisfy BT < ~, in which case M = 0 satisfies 2 M + 1/> 4BT, which yields a small bias, and we can still choose T ~, ro to obtain a small variance for the estimator (7a) (assuming a finite statistical-dependence time ro exists). Thus, both these cases are subsumed by the case of a known finite separable form (11) for the nonstationary correlation. Let us now return to the case of unknown form of nonstationarity. The bias of the estimator (7) is given by bias --a E{R~y( t, "r)}- gxy( t, ~)
= Y. am(z)b,,(t), Iml>M
(24)
and is therefore iclentical to the approximation error in (23). Thus, in order to maintain a small chance of large bias, we require M = 2BT. We shall see in the following that this guarantees that the variance of the estimator (7a) will not decrease as Vol. 15, No, 1, July 1988
38
W.A. Gardner / Nonstationary correlation estimation
T is increased. Alternatively, with M fixed, the variance does decrease linearly in 1/T; however, the chances of a large bias increase with T. We shall see that if T is chosen large enough to yield a small variance (say 20%) for a given value of M, then the chances of a large bias must be high (95%) and the bias can be very large (95%). To prove this is a no-win situation, let us choose a form of nonstationarity that is the most favorable to the estimator (7), namely, that for which only M = 0 would be n e e d e d / f the form of nonstationarity were known. That is, let x(t) be given by
x(t) = c(t)y(t),
(25)
for a close approximation to c(t), M
c(t)-~
~
(26a)
Rxx(t, r) = c(t)c(t - r)Ryy(~'),
(26b)
Ryy(t, ~) = Ryy(Z).
(26c)
To satisfy the persistence condition (1 lc) required for the existence of a consistent estimator, let c(t) be a + 1-valued waveform with zero crossings given by a sample of a Poisson point process. Then the time-average autocorrelation of c(t) is given by (see [4]) -
.
1
O<~t<~T,
(29)
which yields the close approximation (23) with am(r)=amRyy(Z). It can be shown (using Isserlis' formula--see [4]) that the covariance for d,(t) and d,,(z) in (7b) (for jointly Gaussian zero-mean x(t) and y(t)) is given by
cov(,L(r), ,~,.(r)} =-T-17 foTIf b.(t,)b*(t2)
where y(t) is a zero-mean stationary ergodic Gaussian process. Then
Rxy(t, 7") = c(t)Ryy(~'),
ambm(t),
m---M
[Rxx(q, tl -
x
t 2 ) R y y ( t l - 7, t 1 - / 2 )
+ Rxy(tl, tl - t2+ r) X gxy(t2,
t I + 1")] d/I dt2.
t2-
(30) Substitution of (26) into (30) yields
cov{,L(~'), ~m(r)} - T2
b.(tl)b*(t2)c(fi)c(t2) 2
x [R~y(t~ - t2) + R , ~ ( t l - t2 + "r)
× Ryy(tl - t2- r)] dq dt2
s
R ~ ( z ) & llm ~ f ° c ( t ) c ( t - r ) d t
_
1
,
1
T: a - r J~l~l = e -2hIll,
(27)
where B is the average rate of zero crossings. Let B satisfy B~. 1/Zo, where ro is the width of Ryy(r). Then x(t) is not a locally stationary process. Recall that if c(t) were known then we could choose bo(t) = c(t) and (7) would indeed provide a consistent estimator. But with c(t) unknown, we must resort to using an arbitrary orthonormal set. A particularly convenient set is the set of harmonically related sine waves
bin(t) = e i2~m'/T,
(28)
since it is intuitively clear that we need M-~ 2 B T Signal Processing
X C(t+½v)C(t-½v) dt 2 x [Ryy(v) + Ryy(V + ½z)
x Rry(V -½~')] dr.
(31)
Substitution of (28) into (31) yields
cov{,~.(~'), &.(~)} 1 /"
1 fT-~MC(t+½v )
x c(t-½v) e -i2~"-m)t/r dt x [ R E ( v ) + Ry e (v + ½~')
x R~A ~'- ½~')] e'~'c"+m~/" d~,.
(32)
W.A. Gardner/ Nonstationary correlation estimation Since the product waveform c(t +½v)c(t -½v) contains no finite additive sine wave components, then the limit of the inner integral in (32) as T-->oo is zero for n # m and, therefore, limcov{d,(r),~m(r)}=0,
T~cx3
n#m.
(33)
Consequently for large enough T, the variance of R~y(t, r) from (7a) is closely approximated by M
var{/~r(t, z)}~
~
var{am(r)}lb,,(t)l 2.
rn=--M
(34) Using n = rn in (32) yields var{~,~ ('r)}
If T If T+' C( t + v/2)C( t-- V/2) dt
=-T - r - T J~1~1
× [R~e(Vl+ Rev(v + ~ ) x Ryy(V-½r)] e i2~''0~/r dr.
(35)
Since c(t) has bandwidth B, the inner integral in (35), which closely approximates the time-average autocorrelation (27) of c(t) for T~. I/B, has width on the order of 1/B, whereas [R2ye(V)+ Rye(V+ r / 2 ) R y y ( V - r / 2 ) ] has width on the order of the width of Rye(Z), which is to. Therefore, since B ~ 1~to, (35) is closely approximated by
39
It should be emphasized that this result is indeed exemplary of the behavior of the estimator (7). The result is not due to the particular choice of sine waves for the basis functions, or of the binaryvalued waveform for c(t), or of the Gaussian model for y(t). Rather, it is due to the M >i 2BT rule which applies in general for nearly bandlimited functions Rxy(t, r), -oo < t < oo, not contained in (or closely approximable by functions contained in) a fixed, known, finite-dimensional space. To give the proposed estimator (7) every possible chance of being useful, we could propose that there is somehow a way to reduce the number of basis functions {bin} below 2 M + I ~ 4 B T and thereby obtain a small variance and a tolerable bias. But it is easily seen that without knowledge of the form of nonstationarity, this could only be a very lucky accident, and could not be knowable. Specifically, the L2[0, T] norm of the bias (24) for each r is given by Ilbiasll = Rxy --
m Y~
am(z)bm(t)
m=-M
m~
I
(38)
t
M0
Mo
+
~
am(r)bm(t)
m=M+l
var{8,,(r)}-~-~T[R~e(O)+ R2ye(½r)].
(36)
+ Thus, the variance of each component in the estimator (7a) is inversely proportional to T and does, therefore, approach zero as T-->oo. However, substitution of (36) into (34) yields (using Ib,.(t)J2=l for the complex sine waves and 2 2 Ryy(0) = R~y(t, 0)) the close approximation var{/~y(t, r)} 2 R~y( t, O)
-M-I
am(r)bm(t)
IlbiaslL = 2
lam(,)l 2
= 2
(37)
Now, if we require M >12 B T for small bias, then this normalized variance exceeds 4 no matter how large T is.
(39)
where Moa--2BT>M. Since the most favorable case is for (23) with M = Mo to be exact, let us proceed with this case, in which (39) reduces to
m
2 1 2 ~- (- 2 M - + 1) [1 + Ryy(~r) / Rye(0) ]. BT
am(r)bm(t) I ,
~,
m~--M o
(40)
+1
Now, in order to obtain an acceptably small variance, say 20%, var{kx/t, r)} 2 Rxy( t, O)
1 5 Vol. 15, No. I, July t988
40
W.A. Gardner/ Nonstationary correlation estimation
(37) reveals that we require M<~Mo. Thus, the squared norm of the bias (40) contains the energy of approximately ~19(95%) of the orthogonal components of Rxy. There is little hope of this energy being an acceptably small percentage of the total energy of Rxy. This could only be a very lucky accident (5% chance), and could not be knowable if the form of nonstationarity, c(t), is unknown.
4. Time-variant autoregressive model fitting It is well known that the problem of obtaining a least-squares fit of an autoregressive model to a time-series is equivalent to the problem of estimating the correlation function (sequence) and then solving a set of linear equations (the normal equations) specified by the estimated correlation function. Corresponding to the two methods for estimating nonstationary correlation functions based on (7) and (17), there are two methods for fitting time-variant autoregressive models. When the direct method related to (17) is used, the total record of the time-series to be modeled is segmented, and the procedure for time-invariant model fitting is applied in sequence to each segment. These segments are indexed by t' in (17) and the parameter T in (17) is the segment length. When the indirect method related to (7) is used (see [21, 11 and 18]) T represents the total record length rather than the segment length. Based on this difference it has been argued (see, for example, [6]) that because of the longer averaging time T (over which the squared errors to be minimized are averaged), the method based on the expansion (7) should be superior to the direct method based on (17). But for the problem of estimating correlation functions, the fact that the variances of the individual terms in the expansion (7) add essentially counteracts the reduction in the variance of each term due to the use of an increased averaging time, since the number of terms must be increased. In fact, if N non-overlapping segments of length T= To/N are used in (17), and N terms with T = To are used in (7), then the results in Section Signal Processing
3 reveal that the variances of the autocorrelation estimates for both methods are proportional to N~ To. The lack of superiority of the model-fitting method related to the expansion (7) is demonstrated experimentally for speech modeling by Hall, Oppenheim and Willsky [7]. On the other hand, the utility of the method related to (7) for periodically time-variant autoregressive model fitting and the related problem of linear prediction has been reported by Hasselman and Barnett [8], Parzen and Pagano [16], and Jones and Brelsford [9]. This utility is to be expected for reasons analogous to those given in Section 2 for the problem of estimating periodically time-variant correlation functions using (7).
5. Conclusions The only cases in which nonstationary correlation functions can be consistently estimated or at least accurately estimated with the orthogonal series expansion procedure (7) are those for which (i) the form of nonstationarity is known in advance and is sufficiently persistent, (ii) the nonstationarity is periodic or consists of a sum of periodic components with known periods, or (iii) the nonstationarity is slowly fluctuating (locally stationary). It has been proved that in all other cases it is impossible to guarantee that the general form of estimator (7) will provide an accurate approximation to the correlation function. It has also been shown that the chances of it accidentally providing a good approximation must be small, and that it cannot be known whether or not one has obtained a lucky estimate. This disproves the claim in [13] and [14] that the orthogonal series expansion method provides consistent estimators for essentially any form of nonstationarity. The ramifications of the limitations on the estimator (7) to least-squares time-variant autoregressive model fitting have also been discussed, and the suggestion in [6] that model fitting methods based
W.A. Gardner / Nonstationary correlation estimation
on orthogonal series expansions should be superior to those based on data-segmenting for general nonstationarity is disproved.
References
[1] R.A. Boyles and W.A. Gardner, "Cycloergodic properties of discrete-parameter nonstationary stochastic processes", IEEE Trans. Inform. Theory, Vol. IT-29, 1983, pp. 105-114. [2] F. Donati, "Finite-time averaged power spectra", IEEE Trans. Inform. Theory, Vol. IT-17, 1971, pp. 7-16. [3] L.E. Franks, Signal Theory, Prentice-Hall, Englewood Cliffs, N J, 1969. [4] W.A. Gardner, Introduction to Random Processes with Applications to Signals and Systems, Macmillan, New York, 1985. [5] W.A. Gardner, Statistical Spectral Analysis: A Nonprobabilistic Theory, Prentice-Hall, Englewood Cliffs, N J, 1987. [6] Y. Grenier, "Time-dependent ARMA modeling of nonstationary signals", IEEE Trans. Acoust., Speech, and Signal Process., Vol. ASSP-31, 1983, pp. 899-911. [7] M. Hall, A.V. Oppenheim, and A. Willsky, "Time-varying parametric modeling of speech", Proc. IEEE Decision and Control Conf., New Orleans, 1977, pp. 1085-1091. [8] K. Hasselman and T.P. Barnett, "Techniques of linear prediction for systems with periodic statistics", J. Atmospheric Sci., Vol. 38, No. 10, 1981, pp. 2275-2283. [9] R.H. Jones and W.M. Brelsford, "Time series with periodic structure", Biometrika, Vol. 54, 1967, pp. 403-407.
41
[10] J. Kamp6 de Frriet and F.N. Frenkiel, "Correlations and spectra for non-stationary random functions", Math. Comput., Vol. 16, No. 77, 1962, pp. 1-21. [11] L.A. Liporace, "Linear estimation of nonstationary signals", J. Acoust. Soc. Am., Vol. 58, No. 6, 1975, pp. 1288-1295. [12] R.M. Loynes, "On the concept of the spectrum for nonstationary processes", J. Roy. Statist. Soc., Vol. 30, Ser. B, 1968, pp. 1-30. [13] V.Z. Marmarelis, "A single-record estimator for correlation functions of nonstationary random processes", Proc. IEEE, Vol. 69, 1981, pp. 841-842. [14] V.Z. Marmarelis, "Practical estimation of correlation functions of nonstationary Gaussian processes", IEEE Trans. Inform. Theory, Vol. IT-29, 1983, pp. 937-938. [15] E. Parzen, "Spectral analysis of asymptotically stationary b time series", Bull. Inst. Int. Statist., Vol. 39, No. 2, 1962, pp. 87-103. [16] E. Parzen and M. Pagano, "An approach to modeling seasonally stationary time-series", J. Econometrics, Vol. 9, 1979, pp. 137-153. [ 17] M.B. Priestley and H. Tong, "On the analysis of bivariate non-stationary processes", J. Roy. Statist. Soc., Ser. B, Vol. 35, 1973, pp. 153-188. [18] K.C. Sharman and B. Friedlander, "Time-varying modeling of a class of nonstationary signals", Proc. lnternat. Conf. Acoustics, Speech, and Signal Processing, San Diego, 1984, pp. 22.2.1-22.2.4. [19] R.A. Silverman, "Locally stationary random processes", IRE Trans. Inform. Theory, Vol. IT-3, 1957, pp. 182-187. [20] R.L. Snyder, "A partial spectrum approach to the analysis of quasi-stationary time-series", IEEE Trans. Inform. Theory, Vol. IT-13, 1967, pp. 579-587. [21] T. Subba Rao, "The fitting of non-stationary time-series models with time-dependent parameters", J. Roy. Statist. Soc. Ser. B, Vol. 32, 1970, pp. 312-322.
Vol. 15, No. 1, July 1988
31
C O R R E L A T I O N E S T I M A T I O N A N D TIME-SERIES M O D E L I N G FOR N O N S T A T I O N A R Y P R O C E S S E S * William A. G A R D N E R (Member EURASIP) Signal and Image Processing Laboratory, Department of Electrical Engineering and Computer Science, University of California, Davis, CA 95616, U.S.A. Received 18 February 1987 Revised 9 November 1987
Abstract. Various types of statistical signal processing require the estimation of time-variant correlation functions and the fitting of time-variant models for nonstationary processes. In this paper, the situations for which nonstationary probabilistic correlation functions can be accurately estimated from a single sample path of a stochastic process are delineated. These include only (i) known form of nonstationarity, (ii) periodic or almost periodic nonstationarity, and (iii) slowly fluctuating nonstationarity. Known methods of estimation for each of these situations are reviewed within the unifying framework of orthogonal series expansions. The fact that estimators based on orthogonal series expansions for nonstationary correlation functions other than (i)-(iii) cannot be guaranteed to be accurate is established. Ramifications for time-variant autoregressive model fitting are discussed. Zasammenfassung. Verschiedene Aufgabenbereiche der Signalverarbeitung erfordern die Sch~itzung zeitvarianter Korrelationsfunktionen und die zeitvariante Modellierung nichtstation~irer Prozesse. lm folgenden werden die Situationen dargelegt, fiir die nicht-station~ire Korrelationsfunktionen aus einem einzigen Repr~isentanten eines Prozesses gesch~itzt werden k/~nnen. ErfaJ3t werden kSnnen nur F~ille mit (i) bekannter Form der Instationarit§t, (ii) periodischer oder fast-periodischer Instationaritat und (iii) langsam ver/inderlicher Instationarit~it. Die bekannten Sch~itzverfahren fiir jede dieser Gegebenheiten werden in einem einheitlichen Rahmen, dem der Orthogonalreihen-Entwicklung, zusammengestellt. Die Grenzen der Nutzung von Reihenentwicklungen bei der Sch~itzung nicht-station~irer Korrelationsfunktionen werden aufgezeigt. VariationsmSglichkeiten fiir die zeitvariante AR-Modellierung werden diskutiert. Rrsumr. De nombreuses operations de traitement du signal requi~rent l'estimation de fonctions de corrrlation dependant du temps et la construction de modules variant avec le temps pour les processus nonstationnaires. Dans cet article, les cas pour lesquels les fonctions probabilistes de corrrlation nonstationnaires peuvent ~tre estimres de mani~re pr6cise a partir d'une seule trajectoire d'un processus stochastique sont repertorirs. Ces cas incluent seulement: (i) les formes connues de stationnaritr, (ii) les processus p6riodiques ou presque prriodiques, et (iii) les processus nonstationnaires ~ evolution lente. Les m&hodes connues d'estimation pour chacune de ces situations sont reexaminres du point de vue grnrral des reprrsentations en srries orthogonales. Le fait qu'il n'est pas possible de guarantir la prrcision des estimres sous forme de srries orthogonales pour les fonctions de corrrlation nonstationnaires autres que du type (i)-(iii) est 6tabli. Les cons6quences de ce rrsultat pour la construction de modules autorrgressifs drpendant du temps sont examinres. Keywords. Nonstationary correlation, correlation estimation, time-series modeling, nonstationary time-series.
1. Introduction ~ Various types of signal processing require the estimation of time-variant correlation functions and the fitting of time-variant models for nonstationary processes. The most straightforward * This work was supported in part by a grant from ESL, Inc., Sunnyvale. CA.
way to estimate the correlation function of a nonstationary stochastic process using one sample path x(t) over a time interval [0, T] is to use a sliding window average of the lag product l The research reported by W. Music in "Recursive prediction and estimation of periodic nonstationary signals", 1987 IEEE Region 5 Conference, pp. 139-142, duplicates that presented herein, but was mistakenly not attributed to the author of the present paper.
0165-1684/88/$3.50 © 1988, Elsevier Science Publishers B.V. (North-Holland)
32
W..A. Gardner / Nonstationary correlation estimation
x ( t ) x ( t - r ) . When possible, the window must be narrow enough to resolve the nonstationarity but wide enough to provide adequate averaging. An alternative approach proposed by Snyder [20] uses an orthogonal series expansion in t of the lagproduct time-series x ( t ) x ( t - r ) over the entire interval [0, T], with r treated as a parameter. Boyles and Gardner [1] showed that the orthogonal series expansion method yields consistent estimators (mean squared error approaches zero as T-~ oo) for periodic and almost periodic nonstationarity (cyclostationarity) when the complex sinusoids with appropriate frequencies are used as the basis functions, provided that the process exhibits appropriate ergodic properties (cycloergodicity). 2 But the utility of this approach for more general nonstationarity has not been established. This orthogonal series expansion approach can be applied to autoregressive modeling of nonstationary processes by simply substituting the nonstationary correlation estimates for the true nonstationary correlations into the probabilistic normal equations whose solution provides a minimum mean squared error autoregressive model. A related alternative approach uses orthogonal series expansions of the time-variant autoregressive model parameters and then derives empirical normal equations whose solutions provide the least squares estimates of the expansion coefficients. This approach has been proposed and studied by many authors including Jones and Brelsford [9], Subba Rao [21], Liporace [ 11], Hall et al. [7], Parzen and Pagano [16], Hasselman and Barnett [8], Grenier [6], and Sharman and Friedlander [18]. It can be argued that the attainable performance of this alternative approach is tied to that of the former approach. Thus, insight into the general utility of orthogonal series expansions for both nonstationary correlation estimation and nonstationary autoregressive modeling can be gained by studying the former problem which is the more analytically tractable of the two. 2An almostperiodicfunction is a sum of periodic functions with incommensurateperiods. Signal Processing
With the wide-spread applicability of any method for obtaining good estimates of nonstationary correlation functions as motivation, the purposes of this paper are (i) to delineate the types of nonstationarity for which consistent or at least accurate estimates of nonstationary correlation functions can be obtained and to explain how to obtain them, and (ii) to establish that orthogonal series expansions cannot provide estimates of nonstationary correlation functions that are guaranteed to be accurate for any other types of nonstationarity. The remainder of this introduction provides a summary of the current status of the problem of estimation of nonstationary correlation. The conventional method for estimating a correlation function, say Rxy(7. ) a E{x(t)y(t - ~')},
(1)
for a pair of jointly wide-sense stationary stochastic processes x(t) and y(t) using one sample path from each of x(t) and y(t) is to time-average the lag product:
a if ofx ( t ) y ( t - ~') dt.
Rxy(7) =-~
(2)
If x(t) and y(t) possess appropriate ergodic properties (see [4, pp. 150-153] and Section 2), then we have the desired limit A
lim Rxy(r) = Rxy(r)
T~oo
(3)
in some probabilistic sense such as in mean square or with probability one (depending on the particular ergodic property). It has been shown by Parzen [15] that (3) also holds if x(t) and y(t) are only asymptotically jointly wide-sense stationary. In fact, even if x(t) and y(t) are generally nonstationary, but are asymptotically mean stationary in the sense that the time-average Rxy(r) & lim 1 S~oo S
Rxy(t, r) dt
(4)
of the nonstationary correlation function
Rx~,(t, r) a E { x ( t ) y ( t - r)}
(5)
33
W.A. Gardner / Nonstationary correlation estimation
exists and is not identically zero, then (3) holds provided that x(t) and y(t) possess appropriate ergodic properties (see [10], [4, pp. 150-159], and Section 2). More recently, it has been shown by Boyles and Gardner [ 1] that by using an asymptotically orthonormal series expansion the nonstationary correlation function (5) itself (not just its timeaverage) can be consistently estimated (mean square estimation error approaches zero for every t, or almost every t, as averaging time T approaches infinity) provided that the nonstationarity is either cyclostationary or almost cyclostationary (Rxy(t, z) is either periodic or almost periodic--in the mathematical sense---in t) and provided that x(t) and y(t) possess appropriate ergodic properties called cycloergodicity (see also [4, pp. 345-349] and Section 2). It has been claimed by Marmarelis [13, 14] that the orthogonal series method proposed by Synder [20] can provide consistent estimators for more general nonstationary correlation functions, However it is shown here in Section 3 that this is so only if the form of nonstationarity is known in advance.
2. The orthonormal expansion method In order to analyze the consistency of the estimator to be described it must be assumed that the stochastic process models are sufficiently well behaved that the estimator actually exists. This can be accomplished by requiring that the processes x(t) and y(t) be measurable so that almost all sample functions of the lag-product x ( t ) y ( t - T ) are appropriately integrable. Also, in order to prove that the mean squared error of an estimator obtained from an orthonormal expansion approaches zero for every value of t as the integration time T approaches infinity, it must be assumed that R~y(t, r) is a continuous function of t; otherwise only the L 2 norm of the error can be expected to approach zero or the error can be expected to approach zero for every value of t except for a set of zero measure. Of course, for any physical phenomenon for which the correlation estimation
problem makes sense, there will exist an appropriate stochastic process model. Consequently, these technical aspects of the problem are minimized in the treatment presented here since the primary objective is to present arguments and results that are directly useful to the practitioner. Motivated by the classical method of orthonormal expansions of a function, the following method for estimation of a nonstationary correlation function is considered. Let {bm(t)}-MM be a set of orthonormal functions on the interval [0, T], llo r
bm(t)b*(t)dt=
{1,
m=n,
O, m # n.
(6)
The proposed estimator for Rxy(t, r) is given in terms of {b,,(t)}_MM by M
J~xy(/, 7") & ~, t~m(z)bm(t),
(7a)
m=-~ where a,,(T) =-~ ^ ~ llor x(t)y(t-z)b*(t)dt.
(7b)
In the special case for which M = 0 and bo(t) = 1, this method reduces to the conventional method for estimating a stationary correlation function Rxy(t, r)=- Rxy(z), or the time-averaged value,
ao(z)&~l for Rxy(t, ~') dt,
(8)
of a nonstationary correlation function. It is wellknown that the method produces consistent estimators (mean squared error approaches zero as T - * ~ ) in these two cases provided that the stochastic processes x(t) and y(t) exhibit the appropriate ergodic property (see [2] and [4, pp. 150-159]), namely (14) for b,,(t) =- 1. In the special case for which {bm(t)}_MM are sine waves, bm(t)=e i2~m',
ot-,,,=-Otm,
(9)
this method reduces to that studied by Boyles and Gardner [1], where it is shown that it produces consistent estimators provided that: Vol. 15, No. l, July 1988
34
W . A . G a r d n e r / N o n s t a t i o n a r y correlation e s t i m a t i o n
(i) the correlation function R x y ( t , 7") is either periodic ({am}_MM are all integer multiples of a single fundamental frequency) or almost periodic with only M harmonics, M
Rxy(t, 7,)= Y~ am(7.)e i2~ nl,, rn=
M
it can be shown by methods analogous to those used by Boyles and Gardner [1] that the estimator (7) is consistent provided that x(t) and y(t) possess appropriate ergodic properties. This is accomplished by showing that the bias is asymptotically zero, lim E{/~xy(t, ¢¢)}-Rxy(t, 7-)---0,
(12)
a_m(7,)= a*(7.), (lOa) for almost all t and 7-, and that the variance of each one of the coefficient estimators in (7) is asymptotically zero,
where am(~') ~ lim 1 S~oo S
Rxy(t, 7,) e
~ at,
i2 . . . .
(10b)
lim var{~m(7,)}-= 0,
T~oc5
(ii) x(t) and y(t) possess appropriate cycloergodic properties, namely (14) with bin(t) given by (9). In this case, the set {bm(t)}~_Mis not orthonormal unless {a,,}-M M are all integer multiples of a single fundamental frequency and T is an integer multiple of the reciprocal of this fundamental frequency. Nevertheless, the set {bm(t)}_~M is asymptotically (T-->oo) orthonormal and the estimator (7) is indeed consistent as shown by Boyles and Gardner [1] (see also [4, pp. 345-349]). In fact, it is shown by Boyles and Gardner [1] that even if M in (10) is not finite, the limit of the estimators (7) as T ~ o o and M ~ o o is consistent if (14) is satisfied. (For related work on consistent estimation of the Fourier transform of R x y ( t , • ), which is the instantaneous spectrum, see [5].) In the special case for which the correlation function is given by the finite separable form M
Rxy(t, 7.)= ~
a,,(7,)bm(t), O<~t<~,
m=-M
provided that property
x(t) and y(t) possess the ergodic
lfo fo
lim--~
Rz~,(t, u)b*(t)bm(u) dt du
-- 0, where
(13)
(14a)
z~(t) is the centered lag-product process
z~(t) a=x(t)y(t-7,)-E{x(t)y(t-7-)}. (14b) For example, x(t) and y(t) can be jointly Gaussian with correlation functions that asymptotically decay to zero as the lag 7 approaches infinity. (This is demonstrated in the next section.) Of course, the fact that {bm(t)}M_M in (7) and (11) must be one and the same implies that the form of nonstationarity must be known before the estimate (7) can be formed. For example, if
x(t) = bo(t)y(t) and y(t) is stationary, then (1 la) holds with M --- 0, and if bo(t) is persistent in the sense that
(lla) rlim ~ IT
where
for ]bo(t)[ 2 dt = 1,
then (llc) holds, and ( l l b ) becomes
Rxy(t, 7.)b*(t) dt, am(7")~-limTIo r r ~
(llb)
and the functions {bm(t)}_~M are asymptotically orthonormal lim
T~c~
T
lfo
Signal Processing
bm(t)b*(t)dt=
{1,
re=n,
0,
m ~a~ n ,
(llc)
ao(r) = Rry(Z). Furthermore, if y(t) exhibits mean square ergodicity of its autocorrelation (see [4, pp. 150-153]; e.g., y(t) could be Gaussian with no additive random sine wave components), then the estimator (7) with M = 0 is consistent.
W.A. Gardner / Nonstationary correlation estimation
Another example which better illustrates the constraint (11c) is
x(t) = bo(t)Zo(t) + bl(t)zl(t) + b_l(t)z_l(t), where Zo(t), zl(t), z_l(t), and y(t) are zero-mean jointly stationary ergodic Gaussian processes, bo(t)-=l, and bl(t) and b_l(t) are two +1 waveforms whose zero-crossings are specified by two independent samples of a Poisson point process. Then ( l l c ) as well as ( l l a ) and ( l l b ) are satisfied. This illustrates that a wide variety of appropriate {b,,(t)}ffM other than sine waves do indeed exist. The real limitation to the practical applicability of the required form (11) is that the basis functions {b,,(t)}ffM must be known in advance. More generally, if (11a) is only an approximation, but this approximation becomes exact in the limit as M -> oo, then the estimator (7) is consistent in the limit M--> oo, if (14) is satisfied. Or, in other words, for any arbitrarily small positive number e, there exists a positive integer M and a positive number T both sufficiently large to guarantee that the mean squared error of the estimator (7) is less than e for all t e [0, T] and all • such that t - ~"e [0, T]. This is proved for the special case of sinusoidal basis functions (9) by Boyles and Gardner [1]. This suggests that the estimator (7) might be made consistent for any unknown but persistent nonstationarity for which (14) is satisfied by letting M-> oo and using any complete orthonormal set on L2[0, T] (the space of square integrable functions on the interval [0, T]) that is parameterized by T so that (11c) is satisfied. The objectives of the next section are to show that, for general unknown nonstationarity, such an estimator is inconsistent and for finite T and M cannot be guaranteed to be accurate. However, before proceeding, it is explained that for nonstationarity that is neither of known form nor periodic (nor almost periodic), the only situation in which accurate estimates can be guaranteed is that where the nonstationarity is known to be sufficiently slow, typically called local stationarity or quasi-stationarity. It is shown in the next
35
section that in this case (7) is appropriate with only M = O and bo(t)=-1 (although M > 0 can be used as explained by Snyder [20]). Generally speaking, a pair of processes is said to be jointly locally stationary (see [19] and [20]) if the bandwidth B of their nonstationarity and the maximum width ro of their auto- and cross-correlation functions are related by Bro,~ 1.
(15)
The parameter To is a measure of the width of, say, Rxy(t, r) as a function of ~', and B is a measure of the spectral bandwidth of, say, Rxy(t, ~') as a function of t. If the averaging time T in the conventional estimate (2) is chosen to satisfy
T
T~'~o,
(16a, b)
then the mean squared error between /~xy(r) and Rxy(t, r) will be small for 0<~ t ~< T, provided that x(t) and y ( t - s ) are nearly statistically independent for Isl > In order to track the slow nonstationarity, the fixed window [0, T] can be made into a sliding window [t', t ' + T] to obtain the estimate A , allot Rxy(t,'r)=-~
x(t'+t)y(t'+t-z)dt. (17)
The condition (16a) insures a small bias due to time smearing effects, and condition (16b) insures a small variance. (See [20] and, for related work on estimation of the Fourier transform of Rxy(t, • ) for locally stationary processes, see [12, 2, 17 and 5, and references therein].)
3.
Limitations
The fundamental reason why the estimator (7) cannot be made consistent with M ~ oo (and, of course, T - ~ ) for any correlation function Rxy(., .r) that is contained in L 2 [ 0 , T ] for all T < oo is that there is no separable space containing all such functions s. That is, it is not true that there 3 However, there is such a space for all almost periodic functions whose sine wave frequencies form a countable set. Vol. 15, No. 1, July 1988
W.A. Gardner / Nonstationary correlation estimation
36
exists a countable basis {bm(t)}_~oo depending on the parameter T such that for any such function Rxy( ", 7.) and any arbitrarily small positive number e there exists a positive integer M sutticiently large that
Rxy(t, 7.)-
Y~ am(7.)bm(t)
(18)
m=-M
for almost every t c [0, T] for arbitrarily large T, where am(7.) is given by ( l l b ) . What is true is that for a given accuracy e over the interval [0, T], the number M of terms used must grow with T if the basis {bm(t)}_MM is chosen arbitrarily. Consequently, as T is increased in order to reduce the variance of each term (7b) in the estimator (7a), the bias will grow unless M is also increased. But as M is increased, the variance of the sum in (7a) increases. In fact, for any fixed T, as M ~ ce, we obtain from Bessel's equality for a complete orthonormal set {bm(t)}_~oothe limit
lim Rxy(t, 7.)=x(t)y(t--7.),
M~oo
(19)
for almost every t ~ [0, T], which has an unacceptably large variance although the bias is zero. At the opposite extreme, it is shown in this section that for any fixed M, as T ~ , the chances of obtaining a small bias are M out c~, which is zero, although the variance is zero in this case. One approach to quantifying this problem of trading off bias and variance by selection of M and T is described in this section, and is used to illustrate the fact that (7) cannot be guaranteed to provide accurate estimates for general unknown nonstationarity. In order to prove that M must grow with T, we restrict consideration to the subset of all nonstationary correlations for which the nonstationarity is nearly bandlimited to arbitrarily large but finite bandwidth B. We then establish that M must grow linearly with T. (Since this restriction does not impose any limitation to the practical applicability of the result, we do not pursue a more general and necessarily more abstract proof that M must grow at a certain rate with T.) Signal Processing
Although the form of nonstationarity may be unknown, the physics of any practical situation would enable one to set an upper bound on the highest spectral bandwidth B of practical interest. This can be done mathematically by selecting B such that, say, 99% of the energy of every practically feasible function Rxy(., 7.) on the time interval [0, T] in a given practical situation is contained within the spectral band [ - B , B]. The smallest generally useful subset of L2[0, T] for purposes of approximation is defined to be the range space of the operator on L2(-oo, ~ ) that bandlimits to [ - B , B] and then duration limits to [0, T]. Thus, for every Rxy(" , 7") contained in this smallest generally useful subset, there exists a function f~(.) in L2(-oo, co) such that Rxy( ", 7.) is the image of f~(. ) under this operator:
Rxy(t, 7.)=f~_ook(t,u)f~(u)du,
t~[0, r], (20a)
where k(t, given by
u) is the kernel of the operator and is
k(t,u)=g(t)h(t-u),
(20b)
where
O<~t<~ T,
(20c)
otherwise is the duration-limiting part and
h(t - u) -
sin[27.rB(t - u)]
,.(t-u)
(20d)
is the bandlimiting part. Since this operator is compact, its range space is a compact subset of L 2 and is therefore approximately finite-dimensional. That is, there exists a finite dimensional subspace of L2[0, T] such that every function Rxy( ", 7.) contained in this compact subset has, say, 99% of its energy contained in this subspace. Specifically, the orthogonal projection, denoted by Rxy(', 7.), of Rxy(', 7.) onto this subspace satisfies II/~xy(.,
r)-R~y(., 7.)11
IIRxy(', 7.)11
<0.01,
W.A. Gardner/ Nonstationarycorrelationestimation where
Rxv(., z),, a=[T I ~ lRxy(t, z)]2 dt]~
(21)
is the L2[0, T] norm. If this finite dimensional subspace is spanned by the orthonormal basis {b,, (t)} ffM (cf. (6)), then the orthogonal projection is given by M
Rxy(t, z)=
Y~ am(z)bm(t),
(22a)
rn=-M
where
l forRxy(t, z)b*(t) dt.
a,,('r) = ~
(22b)
With the preceding definition in hand, we have the following theorem. Theorem. The smallest generally useful subset of L2[0, T] for approximation of functions with band-
width parameter B has approximate dimensionality D given by D = 4BT. The p r o o f of this theorem follows directly from the classical results (from the theory of functional analysis) on optimum finite-dimensional approximation of compact operators. A practically oriented account is given by Franks [3]. As a consequence of this theorem, the approximation M
E am(T)bm(t) ~Rxy(t, r) m ~ - M
oo
= E
am(z)b,~(t),
rtl = - c x 3
O<-t<~ T, (23) using an arbitrary basis {b,.}_~ for L2[0, T], cannot be expected to be accurate unless 2 M + 1 1>4BT (say M >12BT). In fact, even with M >!2BT (23) will not be an accurate approximation unless an appropriate basis (e.g. each b,.(t) approximately bandlimited to a spectral bandwidth of at most B) is used. To be most favorable to the proposed estimator (7), let us assume that the basis {bin} is
37
appropriate so that M need not be any larger than 2BT. If M is less than 2BT then the fact that a close approximation in (23) could be only a lucky chance insures that an estimate of the form (7) cannot be a close approximation to Rxy(t, z) except as a lucky chance. In fact one could say that the chances of (7) being capable of providing a good estimate are M out of 2BT. This chancy estimation is not of practical value since there is in general no way to know whether you have been lucky or not. Moreover, the chances of being lucky with a fixed M (which is all that can be used in practice) become smaller and smaller as T is increased, as explained in the following. However, before proceeding it is explained that if it is known that R~y(., z) is periodic with period To, then by choosing b,,(t) = exp(i2~rmt/To), (23) will be a close approximation for M=2BTo, regardless of how large T is. Thus, by choosing M = 2BTo to obtain a small bias, and choosing T sufficiently large to render the variance of the estimator (7a) small (assuming the ergodic property (14)), the estimator can be guaranteed to provide accurate estimates. A similar argument applies when Rxy(', r) is almost periodic. Also, when x(t) and y(t) are jointly locally stationary, then we can choose T to satisfy BT < ~, in which case M = 0 satisfies 2 M + 1/> 4BT, which yields a small bias, and we can still choose T ~, ro to obtain a small variance for the estimator (7a) (assuming a finite statistical-dependence time ro exists). Thus, both these cases are subsumed by the case of a known finite separable form (11) for the nonstationary correlation. Let us now return to the case of unknown form of nonstationarity. The bias of the estimator (7) is given by bias --a E{R~y( t, "r)}- gxy( t, ~)
= Y. am(z)b,,(t), Iml>M
(24)
and is therefore iclentical to the approximation error in (23). Thus, in order to maintain a small chance of large bias, we require M = 2BT. We shall see in the following that this guarantees that the variance of the estimator (7a) will not decrease as Vol. 15, No, 1, July 1988
38
W.A. Gardner / Nonstationary correlation estimation
T is increased. Alternatively, with M fixed, the variance does decrease linearly in 1/T; however, the chances of a large bias increase with T. We shall see that if T is chosen large enough to yield a small variance (say 20%) for a given value of M, then the chances of a large bias must be high (95%) and the bias can be very large (95%). To prove this is a no-win situation, let us choose a form of nonstationarity that is the most favorable to the estimator (7), namely, that for which only M = 0 would be n e e d e d / f the form of nonstationarity were known. That is, let x(t) be given by
x(t) = c(t)y(t),
(25)
for a close approximation to c(t), M
c(t)-~
~
(26a)
Rxx(t, r) = c(t)c(t - r)Ryy(~'),
(26b)
Ryy(t, ~) = Ryy(Z).
(26c)
To satisfy the persistence condition (1 lc) required for the existence of a consistent estimator, let c(t) be a + 1-valued waveform with zero crossings given by a sample of a Poisson point process. Then the time-average autocorrelation of c(t) is given by (see [4]) -
.
1
O<~t<~T,
(29)
which yields the close approximation (23) with am(r)=amRyy(Z). It can be shown (using Isserlis' formula--see [4]) that the covariance for d,(t) and d,,(z) in (7b) (for jointly Gaussian zero-mean x(t) and y(t)) is given by
cov(,L(r), ,~,.(r)} =-T-17 foTIf b.(t,)b*(t2)
where y(t) is a zero-mean stationary ergodic Gaussian process. Then
Rxy(t, 7") = c(t)Ryy(~'),
ambm(t),
m---M
[Rxx(q, tl -
x
t 2 ) R y y ( t l - 7, t 1 - / 2 )
+ Rxy(tl, tl - t2+ r) X gxy(t2,
t I + 1")] d/I dt2.
t2-
(30) Substitution of (26) into (30) yields
cov{,L(~'), ~m(r)} - T2
b.(tl)b*(t2)c(fi)c(t2) 2
x [R~y(t~ - t2) + R , ~ ( t l - t2 + "r)
× Ryy(tl - t2- r)] dq dt2
s
R ~ ( z ) & llm ~ f ° c ( t ) c ( t - r ) d t
_
1
,
1
T: a - r J~l~l = e -2hIll,
(27)
where B is the average rate of zero crossings. Let B satisfy B~. 1/Zo, where ro is the width of Ryy(r). Then x(t) is not a locally stationary process. Recall that if c(t) were known then we could choose bo(t) = c(t) and (7) would indeed provide a consistent estimator. But with c(t) unknown, we must resort to using an arbitrary orthonormal set. A particularly convenient set is the set of harmonically related sine waves
bin(t) = e i2~m'/T,
(28)
since it is intuitively clear that we need M-~ 2 B T Signal Processing
X C(t+½v)C(t-½v) dt 2 x [Ryy(v) + Ryy(V + ½z)
x Rry(V -½~')] dr.
(31)
Substitution of (28) into (31) yields
cov{,~.(~'), &.(~)} 1 /"
1 fT-~MC(t+½v )
x c(t-½v) e -i2~"-m)t/r dt x [ R E ( v ) + Ry e (v + ½~')
x R~A ~'- ½~')] e'~'c"+m~/" d~,.
(32)
W.A. Gardner/ Nonstationary correlation estimation Since the product waveform c(t +½v)c(t -½v) contains no finite additive sine wave components, then the limit of the inner integral in (32) as T-->oo is zero for n # m and, therefore, limcov{d,(r),~m(r)}=0,
T~cx3
n#m.
(33)
Consequently for large enough T, the variance of R~y(t, r) from (7a) is closely approximated by M
var{/~r(t, z)}~
~
var{am(r)}lb,,(t)l 2.
rn=--M
(34) Using n = rn in (32) yields var{~,~ ('r)}
If T If T+' C( t + v/2)C( t-- V/2) dt
=-T - r - T J~1~1
× [R~e(Vl+ Rev(v + ~ ) x Ryy(V-½r)] e i2~''0~/r dr.
(35)
Since c(t) has bandwidth B, the inner integral in (35), which closely approximates the time-average autocorrelation (27) of c(t) for T~. I/B, has width on the order of 1/B, whereas [R2ye(V)+ Rye(V+ r / 2 ) R y y ( V - r / 2 ) ] has width on the order of the width of Rye(Z), which is to. Therefore, since B ~ 1~to, (35) is closely approximated by
39
It should be emphasized that this result is indeed exemplary of the behavior of the estimator (7). The result is not due to the particular choice of sine waves for the basis functions, or of the binaryvalued waveform for c(t), or of the Gaussian model for y(t). Rather, it is due to the M >i 2BT rule which applies in general for nearly bandlimited functions Rxy(t, r), -oo < t < oo, not contained in (or closely approximable by functions contained in) a fixed, known, finite-dimensional space. To give the proposed estimator (7) every possible chance of being useful, we could propose that there is somehow a way to reduce the number of basis functions {bin} below 2 M + I ~ 4 B T and thereby obtain a small variance and a tolerable bias. But it is easily seen that without knowledge of the form of nonstationarity, this could only be a very lucky accident, and could not be knowable. Specifically, the L2[0, T] norm of the bias (24) for each r is given by Ilbiasll = Rxy --
m Y~
am(z)bm(t)
m=-M
m~
I
(38)
t
M0
Mo
+
~
am(r)bm(t)
m=M+l
var{8,,(r)}-~-~T[R~e(O)+ R2ye(½r)].
(36)
+ Thus, the variance of each component in the estimator (7a) is inversely proportional to T and does, therefore, approach zero as T-->oo. However, substitution of (36) into (34) yields (using Ib,.(t)J2=l for the complex sine waves and 2 2 Ryy(0) = R~y(t, 0)) the close approximation var{/~y(t, r)} 2 R~y( t, O)
-M-I
am(r)bm(t)
IlbiaslL = 2
lam(,)l 2
= 2
(37)
Now, if we require M >12 B T for small bias, then this normalized variance exceeds 4 no matter how large T is.
(39)
where Moa--2BT>M. Since the most favorable case is for (23) with M = Mo to be exact, let us proceed with this case, in which (39) reduces to
m
2 1 2 ~- (- 2 M - + 1) [1 + Ryy(~r) / Rye(0) ]. BT
am(r)bm(t) I ,
~,
m~--M o
(40)
+1
Now, in order to obtain an acceptably small variance, say 20%, var{kx/t, r)} 2 Rxy( t, O)
1 5 Vol. 15, No. I, July t988
40
W.A. Gardner/ Nonstationary correlation estimation
(37) reveals that we require M<~Mo. Thus, the squared norm of the bias (40) contains the energy of approximately ~19(95%) of the orthogonal components of Rxy. There is little hope of this energy being an acceptably small percentage of the total energy of Rxy. This could only be a very lucky accident (5% chance), and could not be knowable if the form of nonstationarity, c(t), is unknown.
4. Time-variant autoregressive model fitting It is well known that the problem of obtaining a least-squares fit of an autoregressive model to a time-series is equivalent to the problem of estimating the correlation function (sequence) and then solving a set of linear equations (the normal equations) specified by the estimated correlation function. Corresponding to the two methods for estimating nonstationary correlation functions based on (7) and (17), there are two methods for fitting time-variant autoregressive models. When the direct method related to (17) is used, the total record of the time-series to be modeled is segmented, and the procedure for time-invariant model fitting is applied in sequence to each segment. These segments are indexed by t' in (17) and the parameter T in (17) is the segment length. When the indirect method related to (7) is used (see [21, 11 and 18]) T represents the total record length rather than the segment length. Based on this difference it has been argued (see, for example, [6]) that because of the longer averaging time T (over which the squared errors to be minimized are averaged), the method based on the expansion (7) should be superior to the direct method based on (17). But for the problem of estimating correlation functions, the fact that the variances of the individual terms in the expansion (7) add essentially counteracts the reduction in the variance of each term due to the use of an increased averaging time, since the number of terms must be increased. In fact, if N non-overlapping segments of length T= To/N are used in (17), and N terms with T = To are used in (7), then the results in Section Signal Processing
3 reveal that the variances of the autocorrelation estimates for both methods are proportional to N~ To. The lack of superiority of the model-fitting method related to the expansion (7) is demonstrated experimentally for speech modeling by Hall, Oppenheim and Willsky [7]. On the other hand, the utility of the method related to (7) for periodically time-variant autoregressive model fitting and the related problem of linear prediction has been reported by Hasselman and Barnett [8], Parzen and Pagano [16], and Jones and Brelsford [9]. This utility is to be expected for reasons analogous to those given in Section 2 for the problem of estimating periodically time-variant correlation functions using (7).
5. Conclusions The only cases in which nonstationary correlation functions can be consistently estimated or at least accurately estimated with the orthogonal series expansion procedure (7) are those for which (i) the form of nonstationarity is known in advance and is sufficiently persistent, (ii) the nonstationarity is periodic or consists of a sum of periodic components with known periods, or (iii) the nonstationarity is slowly fluctuating (locally stationary). It has been proved that in all other cases it is impossible to guarantee that the general form of estimator (7) will provide an accurate approximation to the correlation function. It has also been shown that the chances of it accidentally providing a good approximation must be small, and that it cannot be known whether or not one has obtained a lucky estimate. This disproves the claim in [13] and [14] that the orthogonal series expansion method provides consistent estimators for essentially any form of nonstationarity. The ramifications of the limitations on the estimator (7) to least-squares time-variant autoregressive model fitting have also been discussed, and the suggestion in [6] that model fitting methods based
W.A. Gardner / Nonstationary correlation estimation
on orthogonal series expansions should be superior to those based on data-segmenting for general nonstationarity is disproved.
References
[1] R.A. Boyles and W.A. Gardner, "Cycloergodic properties of discrete-parameter nonstationary stochastic processes", IEEE Trans. Inform. Theory, Vol. IT-29, 1983, pp. 105-114. [2] F. Donati, "Finite-time averaged power spectra", IEEE Trans. Inform. Theory, Vol. IT-17, 1971, pp. 7-16. [3] L.E. Franks, Signal Theory, Prentice-Hall, Englewood Cliffs, N J, 1969. [4] W.A. Gardner, Introduction to Random Processes with Applications to Signals and Systems, Macmillan, New York, 1985. [5] W.A. Gardner, Statistical Spectral Analysis: A Nonprobabilistic Theory, Prentice-Hall, Englewood Cliffs, N J, 1987. [6] Y. Grenier, "Time-dependent ARMA modeling of nonstationary signals", IEEE Trans. Acoust., Speech, and Signal Process., Vol. ASSP-31, 1983, pp. 899-911. [7] M. Hall, A.V. Oppenheim, and A. Willsky, "Time-varying parametric modeling of speech", Proc. IEEE Decision and Control Conf., New Orleans, 1977, pp. 1085-1091. [8] K. Hasselman and T.P. Barnett, "Techniques of linear prediction for systems with periodic statistics", J. Atmospheric Sci., Vol. 38, No. 10, 1981, pp. 2275-2283. [9] R.H. Jones and W.M. Brelsford, "Time series with periodic structure", Biometrika, Vol. 54, 1967, pp. 403-407.
41
[10] J. Kamp6 de Frriet and F.N. Frenkiel, "Correlations and spectra for non-stationary random functions", Math. Comput., Vol. 16, No. 77, 1962, pp. 1-21. [11] L.A. Liporace, "Linear estimation of nonstationary signals", J. Acoust. Soc. Am., Vol. 58, No. 6, 1975, pp. 1288-1295. [12] R.M. Loynes, "On the concept of the spectrum for nonstationary processes", J. Roy. Statist. Soc., Vol. 30, Ser. B, 1968, pp. 1-30. [13] V.Z. Marmarelis, "A single-record estimator for correlation functions of nonstationary random processes", Proc. IEEE, Vol. 69, 1981, pp. 841-842. [14] V.Z. Marmarelis, "Practical estimation of correlation functions of nonstationary Gaussian processes", IEEE Trans. Inform. Theory, Vol. IT-29, 1983, pp. 937-938. [15] E. Parzen, "Spectral analysis of asymptotically stationary b time series", Bull. Inst. Int. Statist., Vol. 39, No. 2, 1962, pp. 87-103. [16] E. Parzen and M. Pagano, "An approach to modeling seasonally stationary time-series", J. Econometrics, Vol. 9, 1979, pp. 137-153. [ 17] M.B. Priestley and H. Tong, "On the analysis of bivariate non-stationary processes", J. Roy. Statist. Soc., Ser. B, Vol. 35, 1973, pp. 153-188. [18] K.C. Sharman and B. Friedlander, "Time-varying modeling of a class of nonstationary signals", Proc. lnternat. Conf. Acoustics, Speech, and Signal Processing, San Diego, 1984, pp. 22.2.1-22.2.4. [19] R.A. Silverman, "Locally stationary random processes", IRE Trans. Inform. Theory, Vol. IT-3, 1957, pp. 182-187. [20] R.L. Snyder, "A partial spectrum approach to the analysis of quasi-stationary time-series", IEEE Trans. Inform. Theory, Vol. IT-13, 1967, pp. 579-587. [21] T. Subba Rao, "The fitting of non-stationary time-series models with time-dependent parameters", J. Roy. Statist. Soc. Ser. B, Vol. 32, 1970, pp. 312-322.
Vol. 15, No. 1, July 1988