Study with the several resources on Docsity
Earn points by helping other students or get them with a premium plan
Prepare for your exams
Study with the several resources on Docsity
Earn points to download
Earn points by helping other students or get them with a premium plan
This document, authored by mcleod and hipel, discusses the application of exploratory spectral analysis to hydrological time series data. The authors review current methods for estimating univariate spectral density and suggest improvements. They propose that spectral analysis can be seen as a complementary tool to arima model building. The paper introduces a new diagnostic check for arma model adequacy based on the nonparametric spectral density. The authors then perform exploratory spectral analysis on several hydrological time series and suggest some interesting periodicities for further investigation.
Typology: Lecture notes
1 / 17
This page cannot be seen from the preview
Don't miss anything!
Exploratory Spectral Analysis of Hydrological Time Series
A.I. McLeod
Department of Statistical & Actuarial Sciences The University of Western Ontario London, Ontario N6A 5B Canada
K.W. Hipel
Department of Systems Design Engineering University of Waterloo Waterloo, Ontario, Canada N2L 3C
December 1994
2 Exploratory Spectral Analysis
Abstract: Current methods of estimation of the univariate spectral density are reviewed and some improvements are suggested. It is suggested that spectral analysis may perhaps be best thought of as another exploratory data analysis (EDA) tool which complements rather than competes with the popular ARIMA model building approach. A new diag- nostic check for ARMA model adequacy based on the nonparametric spectral density is introduced. Two new algorithms for fast computation of the autoregressive spectral density function are presented. A new style of plotting the spectral density function is suggested. Exploratory spectral analysis of a number of hydrological time series is per- formed and some interesting periodicities are suggested for further investigation. The application of spectral analysis to determine the possible existence of long memory in riverflow time series is discussed with long riverflow, treering and mud varve series. A comparison of the estimated spectral densities suggests the ARMA models fitted previ- ously to these datasets adequately describe the low frequency component. The software and data used in this paper are available by anonymous ftp from fisher.stats.uwo.ca in the directory pub\mhts.
Key words: AIC-Bayes, autoregressive spectral density estimation, diagnostic checks for ARMA models, exploratory data analysis, fast Fourier transform, Hurst coefficient, long-memory time series, periodogram smoothing, riverflow time series, spectral density plots
4 Exploratory Spectral Analysis
Theorem 6, p.210). One necessary condition for this is that
lim n→∞
n
∑^ n
k=
γk = 0. (1)
A simple example of a stationary non-ergodic time series is the symmetrically correlated time series which has TACVF given by
γk = γ 0 , if k = 0,
= γ, if k = 0,
where γ < γ 0. Ergodicity is a theoretical condition which is not normally possible to verify in practice. In the next section, the basic properties of the spectral density function of a time series are reviewed. In §3 of estimating the spectral density function are discussed. In §4, it is shown that the spectral density function provides the best tool to define what is meant by long memory in time series. The application of spectral methods to annual riverflow and other hydrological data is discussed. A new style of plotting the spectral density function is recommended. Spectral analysis is a component in a PC time series package developed by McLeod and Hipel. A student version of this package is available via anonymous ftp as men- tioned in the Abstract.
2 Spectral Analysis Primer
2.1 Spectral Density Function
Spectral analysis can be regarded as the development of a Fourier analysis for sta- tionary time series. Just as in classical Fourier analysis a real function z(t) is repre- sented by a Fourier series, in spectral analysis the autocovariance function of a station- ary time series has a frequency representation in terms of a Fourier transform. This representation was first given by Herglotz (1911) who showed that any positive-definite function, such as the autocovariance function, γk, of a stationary time series can be rep- resented as
γk =
(−π,π]
eiωkdP (ω), (3)
where P (ω) is the spectral distribution function. In the case where a spectral density function exists, we have dP (ω) = p(ω)dω and Herglotz’s equation, eq. (3), can be writ- ten
γk = 2
∫^ π
0
p(ω) cos(ωk)dω. (4)
The function p(ω), −π ≤ ω ≤ π is called the spectral density function and shares many properties of the probability density function. In addition, note that p(ω) is symmet- ric, p(ω) = p(−ω). For mathematical convenience the units of ω are in radians per unit
A.I. McLeod & K.W. Hipel 5
time, which is known as angular or circular frequency. In practice, however, it is more convenient to work in units of cycles per unit time, which is related to ω by the equation ω = 2πf , where f is now in cycles per unit time.
2.2 Periodogram
A natural estimator of the spectral density function given n observations z 1 ,... , zn from a covariance stationary time series, is given by the periodogram,
I(fj ) =
n
∑^ n
t=
zte−^2 πfj^ (t−1)
2 , (5)
where fj = j/n, j = [−(n − 1)/2],... , 0 ,... , [n/2], where [•] denotes the integer part function. Since I(fj ) = I(−fj ) the periodogram is symmetric about 0, and so when the periodogram or spectral density is plotted we only plot the part where fj > 0. When fj = 0, I(0) = n¯z^2 , where ¯z =
∑n t=1 zt/n. This component,^ I(0), is usually very large due to a non-zero mean and so is ignored in the periodogram and spectral plots. In the case where the spectral density function exists, the expected value of I(fj ) can be shown to be approximately equal to p(fj ), where p(f ) is the spectral density function. In fact, in large samples, I(fj ) for j = 1,... , [(n − 1)/2] are statistically in- dependent and exponentially distributed with mean p(fj ). Schuster (1898) developed the periodogram for searching for periodicities in time series. Sometimes I(fj ) is plotted against its period 1/fj but this is not very satisfactory because, since I(fj ) is calculated at equi-spaced frequencies, the low-frequency part is too spread out. A new style of plotting is suggested in §4 which shows the periodicities on the x-axis.
2.3 Frequency Interpretation
Another important property for the interpretation of the periodogram is that it can be shown that I(fj ) is proportional to the square of the multiple correlation between the observed data sequence z 1 ,... , zn and a sinusoid having frequency fj. Specifically, consider the regression,
zt = A 0 + Aj cos(2πfj ) + Bj sin(2πfj ) + et, (6)
where et is the error term. Then the least-squares estimates of A 0 , Aj and Bj are given by
A 0 =
n
∑^ n
t=
zt
Aj =
n
∑^ n
t=
zt cos(2πfj )
Bj =
n
∑^ n
t=
zt sin(2πfj ).
A.I. McLeod & K.W. Hipel 7
where ck denotes the sample autocovariance function given by
ck =
n
∑^ n
=k+
(zt − ¯z)(zt−k − z¯) for k ≥ 0 ,
and for k < 0, ck = c−k.
2.5 Aliasing
There is an upper limit to the highest frequency that can be observed in the time series. This upper limit, which is 0.5 cycles per unit time or π radians per unit time, is called the Nquist frequency. This upper limit arises because of the discrete time nature of our time series. There is no such upper limit in the continuous time case. To see why aliasing occurs, let f˜ denote any frequency in the interval [0, 0 .5] and let f¨ = f˜ + 0.5. Then it is easily shown for all integer t that cos( f¨ ) = cos( f˜ ) and sin( f¨ ) = sin( f˜ ). The frequencies f˜ and f¨ are said to be aliases. Aliased frequencies, such as f¨ , are observa- tionally indistinguishable from frequencies in the range [0, 0 .5].
2.6 Spectral Density and ARMA Model
The ARMA(p, q) model, may be written in operator notation as,
φ(B)zt = θ(B)at,
where, φ(B) = 1 − φ 1 B −... − φpBp, θ(B) = 1 − θ 1 B −... − θq Bq^ , at is white noise with variance σ^2 a and B is the backshift operator on t, is said to be not redundant if and only if φ(B) = 0 and θ(B) = 0 have no common roots. Due to stationarity and invertibility, all roots of the equation φ(B)θ(B) = 0 are assumed to be outside the unit circle. When q = 0 this model is referred to as the autoregression of order p or AR(p). Fitting high order autoregressive models to estimate the spectral density will be discussed in §3. The theoretical spectral density for time series generated by the ARMA(p, q) model is given by p(f ) = σ a^2 |ψ(e^2 πif^ )|^2 ,
where
ψ(B) =
θ(B) φ(B)
3 Spectral Density Estimation Methods
3.1 Periodogram Smoothing
This is the frequency domain approach because in this approach one works with the Fourier transform of the data sequence rather than the original data sequence. It is the most natural approach to estimating p(f ) although the other approach by high- order autoregressive modelling will often produce a more accurate overall estimate in
8 Exploratory Spectral Analysis
the sense of the integrated mean-square error criterion. On the other hand, periodogram smoothing is better suited to revealing bumps and narrow peaks in the spectrum which is often the principal goal of exploratory spectral analysis in the first place. The Discrete Fourier Transform (DFT) of the sequence z 1 ,... , zn, is defined as
Zj =
∑^ n
t=
zte−^2 πitfj^ , j = −[(n − 1)/2],... , 0 ,... , [n/2],
where [•] denotes the integer part and fj = j/n. The frequencies fj are referred to as the Fourier frequencies. The functions cos(2πfj t) and sin(2πfj t) are orthogonal with respect to the usual inner product when evaluated at the Fourier frequencies. The periodogram is given by
I(fj ) =
n
|Zj |^2.
It follows from the orthogonality mentioned above that
[ ∑n/2]
j=−[(n−1)/2]
I(fj ) =
∑^ n
t=
z^2 t.
Thus I(fj ) can be interpreted as analysis of variance of the data. I(fj ) shows the amount of variation due to each frequency component. Since the periodogram is sym- metric about zero, only positive frequencies need be considered. The periodogram smoothing approach to the estimation of p(f ) is based on the fol- lowing two large-sample results:
< I(fj ) >≈ p(fj ),
and cov(I(fj ), I(fk)) ≈ p^2 (fj ), whenj = k,
≈ 0 , whenj = k.
From the above two equations we see that although I(fj ) is an unbiased estimator of p(fj ), it is not consistent. If p(f ) is assumed to be a smooth function of f , an estimator with smaller mean-square error can be obtained by averaging values of the periodogram. The periodogram smoother may be written
pˆ(fj ) =
∑^ i=q
i=−q
wiI(fj+i),
where q is the half-length of the smoother and the weights wi satisfy the following condi- tions: (i) wi ≥ 0,
10 Exploratory Spectral Analysis
This is a popular method since it seems to produce an estimate of the spectral density which often has lower integrated mean square error than the periodogram smoothing approach. Akaike (1979) addressed the problem of choosing the model order p by suggesting that the weighted average of all spectral functions of all autoregressions of orders p = 0 , 1 ,... , K, where K is some upper limit, should be used. Each spectral density in the average is weighted according to the quasi-likelihood or Bayes posterior given by
e−^
(^12) AIC .
This weighted average estimate is itself equivalent to a special AR(K) model with co- efficients determined by the weighting. This model is referred to as the AR-AIC-Bayes filter. Typically, the AR-AIC-Bayes filter produces a spectral density which shows more features similar to the estimate produced by periodogram smoothing than simply fitting a fixed order AR model. As recommended by Percival and Walden (1993, p.414–417) the AR model parame- ters are estimated using the Burg algorithm rather than the standard Yule-Walker esti- mates. The reason for this is that the standard Yule-Walker estimates are now known to be severely biased in some situations.
3.3 Algorithms For Fast Autoregressive Spectral Density Computation
The spectral density function is normally interpreted by plotting it over the inter- val (0, 12 ). This requires its evaluation at a large number of values. We provide two new faster algorithms for this computation. Fast AR spectrum evaluation is of interest when p is large. In applications, large values of p are not uncommon particularly with the AR- AIC-Bayes filter and also in other applications such as in Percival and Walden (1993, p.522). The two algorithms given can easily be extended to the evaluation of the spec- tral density function of ARMA models. The direct method evaluates the spectral density function, f (λ) given the param- eters φ 1 ,... , φp and σ a^2 directly and in general requires O(np) floating point operations (flops) in addition to O(np) complex exponential evaluations. This method may be improved by avoiding the calculation of the large num- ber of complex exponentials. First, note that the evaluation of the complex expo- nential is equivalent to one sine and cosine function evaluation. So if f (λ) is evalu- ated at equally spaced values throughout [0, 12 ] then the necessary trigonometric func- tions may be evaluated recursively using the sum of angles formulae for sine and co- sine functions. This technique has been previously used in the evaluation of the dis- crete Fourier transform at a particular frequency, see Robinson (1967, pp. 64–65). Given the autoregressive parameters φ 1 ,... , φp and innovation variance σ a^2 , the re- cursive algorithm for spectral density computation may be summarized as follows: Step 1: Initializations. Select M the number of equi-spaced frequencies. Typically, M ← 256 is adequate. Set k ← 1, sβ ← 0, cβ ← 1, sm ← sin(π/(m − 1)) and cm ← cos(π/(m − 1)). Step 2: Compute p(k/2(m − 1)) (a) Set j ← 1, sα ← sβ , cα ← cβ , s ← sα c ← cα A ← 0 and B ← 1.
A.I. McLeod & K.W. Hipel 11
(b) A ← A − sφj and B ← A − cφj (c) j ← j + 1 (d) If j > p then p(k/2(m − 1)) ← σ a^2 /(A^2 + B^2 ) and go to Step 3. Otherwise if j ≤ p then t ← c, c ← cαt − sαs, s ← sαt + cαs, and return to (b) Step 3: Increment for next k Set k ← k + 1. Terminate if k > m. Otherwise if k ≤ m then t ← cβ , cβ ← cmt − smsβ , sβ ← smt + cmsβ , and return to Step 2 In many situations an even faster method is based on the fast Fourier transform (FFT). Given a sequence {ψk}, k = 0,... , M the discrete Fourier transform is a se- quence {Ψ}, = 0,... , M given by
k=
ψe
i 2 πk M (^).
When M = 2q^ for some positive integer q, the discrete Fourier transform may be eval- uated using an algorithm known as the FFT which requires only O(M log 2 (M )) flops when M is a power of 2. The FFT may be applied to evaluate f (λ) by setting ψ 0 = −1, ψk = φk, k = 1 ,... , p and ψk = 0 , k > p. To evaluate f (λ) at N = 2 r equi-spaced points on [0, 12 ], set M = 2 N and apply the FFT. Then f ( /(2N )) = σ^2 a/(Ψ( ) Ψ(¯ )), = 0,... , N , where Ψ(¯ ) denotes the complex conjugate of Ψ( ). In many situations, N = 256 is adequate. The FFT method is the most expedient and convenient method to use with higher level programming languages such as Mathematica in which it is best to vectorize the calculations. In such languages the FFT is provided as a system level function. For ex- ample, in Mathematica (Wolfram, 1991) the following function evaluates the spectral density with σ a^2 = 1 at 256 equi-spaced points:
Arspec[phi_] := 1./(512Re[#1Conjugate[#1]] & )[ Take[Fourier[Join[{1}, -phi, Table[0, {511 - Length[phi]}]]], 256]]
In the above Mathematica function the parameter phi is the vector (φ 1 ,... , φp). The recursive algorithm may be preferable to the FFT method in certain situations where the number of frequencies at which the spectral density is to be evaluated is not a power of 2. This is the case, for example, when spectral methods are used to estimate the parameters of ARMA models (Hannan, §VI.5). The three algorithms outlined above were programmed in Fortran. The FFT method used the algorithm of Monro (1976). Timings for N = 256 equi-spaced eval- uations of f (λ) on [0, 12 ] on a 286/7 PC are shown in Table I. For faster computers, the times will be reduced. Counting the number of flops the recursive method requires O(N p) while the FFT method requires O(N log 2 (N )). This suggests as a rough esti- mate, for p > 8 it may be expected that the FFT will be faster. In fact, as shown in Table I, the FFT method is faster for p ≥ 10.
TABLE I
A.I. McLeod & K.W. Hipel 13
4.1 Exploratory Spectral Analysis of Hydrological Time Series
Hipel and McLeod (1978) previously analyzed some of the longest available station- ary hydrological time series to see if there was a long-memory effect which could not be accounted for by short-memory ARMA models. Some of the data from that study is described in Table II. These datasets include the longest available annual unregulated riverflows as well as long time series on water level, mud varves and treerings. If the long-memory hypothesis is correct then one would expect to see high power in the lowest frequencies and the longer the series length, the more pronounced this ef- fect would be expected to become. It is of interest to compare the two methods of ex- ploratory spectral estimation described in §3. With periodogram smoothing it is ex- pected that the estimates of the spectral density will be more variable but perhaps bet- ter able to resolve the peak at the zero frequency. In the autoregressive approach we will set K quite large so that a good approximation is obtained for any possible long memory effect. The spectral plots are shown in Figures 1–20. Panel (a) in each Figure refers to the periodogram smooth and Panel (b) to the AR-AIC-Bayes Filter. The plotting style for the spectral density shows the periods instead of frequencies on the x-axis. But the scaling on the x-axis is equi-spaced on the frequency scale. This scale desirable because every interval on the frequency scale corresponds to the same number of periodogram estimates. All that changes is just the particular labels used on the axis. These labels are more immediately informative since in most applications most practioners think in terms of the period not of the frequency. The y-axis shows the log to the base 2 of the spectral density function estimate. Cleveland (1994, p.122) pointed out that it is often preferable in statistical graphics to use a log to the base 2 transformation for one of the axes since it is much easier to interpret. When log to the base 2 is used a change of one unit means a doubling of the value in the orginal domain. Similarly, an increase of 0.5 on a log to the base 2 scale, implies a 40% increase. On the other hand fractional powers of ten are much harder to interpret. In the bottom left-hand corner of each plot of the periodogram smooth the shape of the smoother is inset. Only the right half of the bell shaped smoother is shown. The ordinate scale is just chosen for convenience and is not important. The width of the bell shows the number of frequencies involved and its shape indicates the relative weights, wi, in the smooth. In general the two spectral estimation methods, periodogram smoothing and the AR-AIC-Bayes Filter provide quite similar estimates of the spectral density. In general, both methods show that their is high power near the zero frequency which is consistent with the long-memory hypothesis. Some interesting differences between these methods occur for mstouis, neumunas, bigcone, naramata, and navajo. In each of these cases, the periodogram smooth suggests a peak close to zero but not exactly at zero whereas the AR-AIC-Bayes Filter has the peak at zero. In some cases the periodogram smooth suggests interesting possible periodicities in the time series such as in the case of the gota or the ninemile.
14 Exploratory Spectral Analysis
Hydrological Time Series.
Code Name Description n
danube Danube River, Orshava, Romania, 1837–1957, annual flow 120 gota Gota River, Sjotorp-Vanersburg, Sweden, 1807–1957, annual flow 150 mstouis Mississippi River, St. Louis, Missouri, 1861–1957, annual flow 96 neumunas Neumunas River, Smalininkai, Russia, 1811–1943, annual flow 132 ogden St. Lawrence River, Ogdensburg, N.Y., 1860–1957, annual flow 96 rhine Rhine River, Basel, Switzerland, 1807–1957, annual flow 150 minimum Nile River, Rhoda, Egypt, 622–1469, annual minimum level 848 espanola Espanola, Ontario, -471– -820, mud varve 350 bigcone Big cone spruce, Southern California, 1458–1966, treering width 509 bryce Ponderosa pine, Bryce, Utah, 1340–1964, treering width 625 dell Limber pine, Dell, Montana, 1311–1965, treering width 655 eaglecol Douglas fir, Eagle Colorado, 1107–1964, treering width 858 exshaw Douglas fir, Exshaw, Alberta, 1460–1965, treering width 506 lakeview Ponderosa pine, Lakeview, Oregon, 1421–1964, treering width 544 naramata Ponderosa pine, Naramata, B.C., 1451–1965, treering width 515 navajo Douglas fir, Belatakin, Arizona, 1263–1962, treering width 700 ninemile Douglas fir, Nine Mile Canyon, Utah, 1194–1964, treering width 771 snake Douglas fir, Snake River Basin, 1281–1950, treering width 669 tioga Jeffrey pine, Tioga Pass, California, 1304–1964, treering width 661 whitemtn Bristlecone pine, California, 800–1963, treering width 1164
4.2 Comparison with the ARMA Models
For comparison, Panel (c) in each Figure shows the fitted spectral density for the ARMA models estimated by Hipel and McLeod (1978). The specific ARMA model is indicated on the plot. These ARMA models were fitted using the approximate likelihood algorithm of McLeod (1977). In many cases, such as for minimum and espanola, the best fitting ARMA has even higher spectral mass near zero than either of the spectral methods. Only in the case of the ninemile does the spectral density function from the fitted ARMA model seem to be drastically out of line. In this case, after refitting, it was found that an ARMA(9,1) model with φ 2 =... = φ 8 = 0 provided a better fit. The new diagnostic plot is shown in Figure 21. The period of about ten years for the ninemile series is very close to the annual sunspot numbers. Using a pre-whitening analysis, it is found that there is a rather large negative correlation, -0.121 (± 0 .06 sd) between treering width and sunspot activity in the previous year. It is known that high sunspot activity is strongly related to the incidence of melanoma in human males in Conneticut two years later (Andrews and
16 Exploratory Spectral Analysis
References
Akaike, H. (1979), A Bayesian extension of the minimum AIC procedure of autore- gressive model fitting. Biometrika 66(2), 237–242.
Andrews, D.F. and Herzberg, A.M. (1984), Data: A Collection of Problems from Many Fields for the Student and Research Worker , New York: Springer-Verlag.
Barnard, G.A. (1956), Discussion on Methods of using long term storage in reser- voirs, Proc. Inst. Civil Eng. 1, 552–553.
Beran, J. (1992), Statistical methods for data with long-range dependence, Statisti- cal Science 7(4) 404–427.
Bloomfield, P. (1976), Fourier Analysis of Time Series: An Introduction, New York: Wiley.
Box, G.E.P. and Jenkins, G.M. (1975) Time Series Analysis: Forecasting and Con- trol, (2nd edn). San Francisco: Holden-Day.
Brillinger, D.R. (1981), Time Series: Data Analysis and Theory, Expanded Edition, San Francisco: Holden-Day.
Brillinger, D.R. and Krishnaiah, P.R. (1983), Time Series in the Frequency Domain, Amsterdam: North Holland.
Cleveland, W.S. (1994), The Elements of Graphing Data, Revised Edition, Summit: Hobart Press.
Granger, C.W.J. and Joyeux, R. (1980), An introduction to long memory time se- ries models and fractional differencing. Journal of Time Series Analysis 1(1), 15–29.
Hamming, R.W. (1977), Digital Filters, New Jersey: Prentice-Hall.
Hannan, E.J. (1970), Multiple Time Series, New York: Wiley.
Herglotz, G. (1911), Uber potenzreihen mit positivem reellem teil im einheitskreis.¨ Sitzgxber. Sachs Akad. Wiss. 63, 501–511.
Hipel, K.W. & McLeod, A.I. (1994), Time Series Modelling of Water Resources and Environmental Systems, Amsterdam: Elsevier.
Hosking, J.R.M. (1981), Fractional differencing, Biometrika, 68(1), 165–176.
Mandelbrot, B.B. and Van Ness, J.W. (1968), Fractional Brownian motion, frac- tional noises and applications. SIAM Review 10(4), 422–437.
Mandelbrot, B.B.and Wallis, J.R. (1969), Some long-run properties of geophysical records. Water Resources Research 5(2), 321–340.
McLeod, A.I. (1977), Improved Box-Jenkins estimators. Biometrika 64(3), pp.531–
A.I. McLeod & K.W. Hipel 17
Monro, D.M. (1976) Algorithm AS 97. Real discrete fast Fourier transform. Applied Statistics 25(2), 173–180.
Percival, D.B. and Walden, A.T. (1993) Spectral Analysis for Physical Applications, Cambridge University Press.
Priestley, M.B. (1981), Spectral Analysis and Time Series, New York: Academic Press.
Schuster, A. (1898), On the investigation of hidden periodicities with application to a supposed 26 day period of meteorological phenomena. Terr. Magn. 3, 13–41.
Robinson, E.A. (1967). Multichannel Time Series Analysis with Digital Computer Programs, San Francisco: Holden-Day.
Tukey, J. W. (1978), Comment on a Paper by H.L. Gray et al. in The Collected Works of John W. Tukey, Volume II, Time Series: 1965–1984. Ed. D.R. Brillinger. Monterey: Wadsworth.
Wolfram, S. (1991) Mathematica: A System for Doing Mathematics by Computer, (2nd edn). Redwood City: Addison-Wesley.
