Notes for Tuesday, April 14, 1998

Cyclic Behavior of Time Series: the Spectrum and the Periodogram

Many times series show cyclic fluctuations (hours, days, weeks, months, seasons, years), but the underlying sine and cosine functions to use to model this are usually not obvious.

One way to approach this problem is to look at the discrete Fourier transform (DFT) of the time series. The sequence of squared Fourier coefficients in the DFT is called the spectrum of the time series. The periodogram is an estimate of the spectrum, based on the DFT of the observed time series, that can help us to see what the important cycles in a time series are.

This sort of analysis is called “spectral analysis”, or “looking at the time series in the frequency domain.” We will first look at the spectrum and periodogram and spectrum from the point of view of Fourier analysis, and then we will look at an interesting connection with linear regression. Finally we will look at some examples.

Short Digression: Fourier Analysis of f(t), t ∈ [0, 1]

A standard theorem from differential equations, advanced calculus, or even a sophomore/junior linear algebra class, says that any reasonable real function f(t) with domain [0, 1] (or periodic function with period 1) has a Fourier series representation,

f(t) = ∑_{k=1}^{∞} a_{k} cos(2πkt) + b_{k} sin(2πkt) = ∑_{k=−∞}^{∞} c_{k}e^{2πikt},

where a_{k} and b_{k} are real coefficients, e^{iθ} = cos(θ) + i sin(θ), c_{k} = (a_{k} − ib_{k})/2, and c_{−k} = \overline{c_{k}} = (a_{k} + ib_{k})/2. The coefficients are obtained by

a_{k} = 2 ∫_{0}^{1} f(t) cos(2πkt)dt and b_{k} = 2 ∫_{0}^{1} f(t) sin(2πkt)dt  (∗)

or equivalently c_{k} = ∫_{0}^{1} f(t) exp(2πikt)dt. The set of coefficients {a_{k}, b_{k} : k = 1, 2, …} or equivalently {c_{k} : k = 1, 2, …} is called the Discrete Fourier Transform (DFT) of f(t).

Recalling that for any complex number z = a + ib, z\overline{z} = (a + bi)(a − bi) = a^{2} + b^{2} = |z|^{2}, we also know that

∫_{0}^{1} f^{2}(t)dt = ∫_{0}^{1} f(t)\overline{f(t)}dt = ∫_{0}^{1} f(t) ∑_{k=−∞}^{∞} c_{k}e^{−2πikt} dt = ∑_{k=−∞}^{∞} c_{k}\overline{c_{k}} = ∑_{k=−∞}^{∞} |c_{k}|^{2} = 2 ∑_{k=1}^{∞}(a_{k}^{2}+b_{k}^{2})/4

this is called Parseval’s identity. The function

S(k) = |c_{k}|^{2} = (a_{k}^{2} + b_{k}^{2})/4

is called the spectrum of f(t); a particular estimate of S(k) from time series data, described below, is called a periodogram.

The point: If S(k) ≈ 0 for all k beyond some K, we might as well truncate the Fourier series for f(t) to just K terms.

Fourier Analysis of a Time Series as Linear Regression

This is not quite the way spectral analysis of time series is usually developed, but it will give you the right intuition for the problem.

Consider a very simple linear regression model for a time series y_{t} exhibiting cyclic variation with known period p, and overall mean ∑_{t=1}^{n} y_{t}/n ≈ 0:

y_{t} = α cos(ωt) + β sin(ωt) + z_{t},

where z_{t} is a white noise sequence and ω = 2π/p is the known frequency, i.e. [after some trig identities; try fooling around with e^{i(a+b)}] we are fitting a phase-shifted sinusoid

\sqrt{α^{2} + β^{2}} sin(ωt + t_{0})

to y_{t}. This may be rewritten in matrix form as

[ y_{1} ; y_{2} ; ⋮ ; y_{n} ] = [ \begin{matrix}cos(ω1) & sin(ω1) \\ cos(ω2) & sin(ω2) \\ ⋮ & ⋮ \\ cos(ωn) & sin(ωn) \end{matrix} ] [ \begin{matrix}α \\ β \end{matrix} ] + [ z_{1} ; z_{2} ; ⋮ ; z_{n} ] , or \underline{y} = X\underline{β} + \underline{z}

We may use the usual matrix formula for the regression coefficients ˆ \underline{β} = (X′X)^{−1}X′\underline{y} together with some trig integrals from Calc II, to show that

α̂ ≈ \frac{2}{n} ∑_{t=1}^{n} y_{t} cos(ωt)

β̂ ≈ \frac{2}{n} ∑_{t=1}^{n} y_{t} sin(ωt)

Clearly, these are discrete approximations to the integrals in (∗) defining a one-term DFT for f(t) = y_{t/p}.

Moreover, the regression sum of squares is

∑_{t=1}^{n} ŷ_{t}^{2} = ∑_{t=1}^{n}[α̂ cos(ωt) + β̂ sin(ωt)]^{2} ≈ \frac{n}{2}(α̂^{2} + β̂^{2})

(a version of Parseval’s identity), so that R^{2} = ∑_{t} ŷ_{t}^{2}/ ∑_{t} y_{t}^{2} is determined by the squared amplitude of the fitted sinusoid.

This is a quite simplistic model. More realistic would be if we not only did not know the amplitude but also didn’t know the frequency of the series.

So consider a model of the form

y_{t} = ∑_{k=1}^{K} {α_{k} cos(ω_{k}t) + β_{k} sin(ω_{k}t)} + z_{t};

we are now fitting K shifted sinusoids to y_{t}. This is again an easy model to estimate, provided we pick the right set of frequencies. If we take ω_{k} = 2πk/n, with k < n/2 then some nice things happen:

• For starters, this is the only set of frequencies that we really can estimate: smaller periods would wiggle more often than the data, longer periods would not see an entire cycle in the data!

• By fooling around with x = e^{2πikt/n} in the identity (1 − x^{n}) = (1 − x)(1 + x + x^{2} + ⋯ + x^{n−1}), one can show that e^{2πikt/n}, k = 1, 2, … represents a set of orthogonal functions. That is the correlation between any pair of these functions (measured at t = 1, 2, …, n) is zero.

• The nicest thing that happens (because of the orthogonality) is that the estimates of α̂ and β̂ are the same as before, namely,

α̂_{k} = \frac{2}{n} \{∑_{t=1}^{n} y_{t} cos(ω_{k}t)\}

β̂_{k} = \frac{2}{n} \{∑_{t=1}^{n} y_{t} sin(ω_{k}t)\}

• Some further manipulations show that

∑_{t=1}^{n}[y_{t} cos(ω_{k}t) + y_{t} sin(ω_{k}t)]^{2} = \frac{n}{2}(α̂_{k}^{2} + β̂_{k}^{2})

so that α̂_{k}^{2} +β̂_{k}^{2} again gives the sums of squares associated with each of the K fitted sinusoids.

• Finally, (via some math, or common sense), if all the frequencies are fit, then the regression model fits exactly.

The Periodogram

Now we define the periodogram function I(ω) for 0 ≤ ω ≤ π, by

I(ω) = \frac{1}{n} ∑_{t=1}^{n}[y_{t} cos(ωt) + y_{t} sin(ωt)]^{2}.

Of course I(2πk/n) = (α̂_{k}^{2} + β̂_{k}^{2})/2 is an estimate of (twice) the spectrum S(k) = (α̂_{k}^{2} + β̂_{k}^{2})/4 of the time series, and is proportional to the squared amplitude of the sinusoid

\sqrt{α_{k}^{2} + β_{k}^{2}} sin(ω_{k}t + t_{k})

fitted at frequency ω_{k} = 2πk/n, and represents that power (amount of energy) that can be estimated at each of these fundamental frequencies.

Notice that we can “parition” the total sum of squares for the data y_{t} as

∑_{t=1}^{n} y_{t}^{2} = nI(0) + 2n ∑_{k=1}^{K_{max}} I(2πk/n)

where I(0) corresponds to the intercept ∑_{t=1}^{n} y_{t} if it is nonzero, and K_{max} = ⌊n/2⌋ is the largest integer less than or equal to n/2. On the one hand, this is just Parseval’s identity again. On the other, it suggests the following ANOVA table, for 1 ≤ K ≤ K_{max}, when ∑ y_{t} = 0:

k  SS  Periodogram
1  SS_{1} = ∑_{t=1}^{n}[y_{t} cos(ω_{1}t) + y_{t} sin(ω_{1}t)]^{2} = 2n I(2π/n)
2  SS_{2} = ∑_{t=1}^{n}[y_{t} cos(ω_{2}t) + y_{t} sin(ω_{2}t)]^{2} = 2n I(2π2/n)
⋮  ⋮  ⋮
K  SS_{K} = ∑_{t=1}^{n}[y_{t} cos(ω_{K}t) + y_{t} sin(ω_{K}t)]^{2} = 2n I(2πK/n)

Regression SSR =  ∑_{t=1}^{n} ŷ_{t}^{2} = ∑_{k=1}^{K} SS_{k}  = 2n ∑_{k=1}^{K} I(2πk/n)

Error  SSE = ∑_{1}^{n}(y_{t} − ŷ_{t})^{2} = ∑_{k=K+1}^{n/2} SS_{k} = 2n ∑_{k=K+1}^{n/2} I(2πk/n)

Total  SST = ∑_{1}^{n} y_{t}^{2} = ∑_{1}^{n} ŷ_{t}^{2} + ∑_{1}^{n}(y_{t} − ŷ_{t})^{2} = 2n ∑_{k=1}^{K} I(2πk/n)

If you go back to your old linear regression notes, and examine the matrix form of regression \underline{y} = X\underline{β} + \underline{ε} when the columns of X are orthogonal, you will see that this is the same ANOVA table: each SS is the dot product of \underline{y} with the corresponding column(s) of X, and (because of orthogonality) it doesn’t matter what order the independent variables are entered in the regression.

Some Examples

Male and Female Respiratory Deaths in the UK. First, we re-plot the male and female deaths in the UK from various respiratory illnesses.

In the male sequence, there appears to be a clear 12-month cycle, and a hint of a cycle at about five or six months. Not much else is evident from these plots.

The DFT of a discrete sequence y_{t} can be conveniently computed in SPLUS with the so-called “fast Fourier transform”, fft(). Thus, a really raw periodogram may be calculated as follows (the Mod function computes the modulus of a complex number):

402 > par(mfrow=c(4,1))
402 > length(uk.male)
[1] 72
402 > plot((1:35)/6,10/log(10)*log(Mod(fft(uk.male)))[2:36],type="l")
402 > title("Really raw periodogram")

The values of the periodogram I(ω_{k}) have been plotted on a decibel scale, 10 · log_{10} S(ω_{k}), since that is how the periodogram is usually plotted.

The overall mean is not zero here, so we should take it out before we compute the periodogram. Often, not only a mean but a linear trend is taken out of the data before calculating the periodogram, for example:

402 > uk.resids _residuals(lm(resp ~ time, data.frame(time=1:72,resp=c(uk.male))))
402 > plot((1:35)/6,10/log(10)*log(Mod(fft(uk.resids)))[2:36],type="l")
402 > title("Really raw periodogram, detrended")

The function spectrum in SPLUS also computes a periodogram from the detrended time series, but in addition it smooths the observations at the beginning and end of the series; and allows us to run one or more moving-average smooths over the final periodogram.

402 > spectrum(uk.male)
402 > spectrum(uk.male,spans=c(3,3))
402 > # smooth with running means of length 3, twice

Graphs of these various periodogram estimates follow. By convention, all results are plotted on a decibel (10 log_{10}) scale.

Notice that the main power corresponds to a 12 months cycle, but there is also significant power at some other frequencies (a half- or one-third-year cycle and a two-year cycle). Remember, this is like looking at an ANOVA table to decide which terms to keep in the regression.

Raw and smoothed periodograms of the uk.female series are essentially indistinguishable from these.

Lutenizing hormone. We will look at the late1 series from these data.

402 > ts.plot(late1)
402 > title("Lutenizing Hormone, Late1 Series")
402 > length(late1)
[1] 48
402 > l1.resids _residuals(lm(resp~time,data.frame(resp=c(late1),time=1:48)))
402 > plot((1:23)/48,10/log(10)*log(Mod(fft(l1.resids)))[2:24],type="l")
402 > title("Really raw periodogram, detrended")
402 > spectrum(late1)
402 > spectrum(late1,span=3) # follow with a single running mean of length 3

The plots below suggest there is a lot of power at three wavelengths, 8 = .17*48, 14 = .29*48 and 19 = .40*48, although there appears to be “junk” at shorter wavelengths as well.

Serial Dependence in Time Series

The second popular way to analyze time series data is “state space analysis” or “looking at the time series in the space domain”. This method concentrates on the serial dependence rather than the cyclic behavior of the series. The most widely used models are ARMA and ARIMA processes. We introduce the basic AR and MA processes here.

Autocovariance and autocorrelation function

Recall that the autocovariance function of a stationary random function {Y (t)} is

γ(k) = Cov{Y (t), Y (t − k)}

This does not not depend on t because of the the assumption of stationarity. Since γ(0) is the variance, the autocorrelation function is

ρ(k) = γ(k)/γ(0)

Basic time series models: AR and MA

The two basic building-blocks for the “noise” part of a time series model are the autoregressive (AR) and moving average (MA) processes. MA processes are easier to describe first.

• A \underline{moving average} process of order q, denotes MA(q), has the form

y_{t} = μ + β_{0}z_{t} + β_{1}z_{t−1} + … + β_{q}z_{t−q}

or, in terms of the backshift operator Bz_{t} = z_{t−1},

y_{t} − μ = (β_{0} + β_{1}B + … β_{q}B^{q})z_{t}
≡ Θ(B)z_{t},

where z_{t} is a mean-zero white noise process. It is easy to show that y_{t} is always stationary, and

γ(k) = \{ \begin{matrix}∑_{t=0}^{q−|k|} β_{t}β_{t+|k|} & |k| ≤ q \\ 0 & |k| > q \end{matrix}

Here are plots of an MA(1) and an MA(2) process.

The autocorrelation plot for a MA(q) process drops to zero (except for noise) after lag q, which make the acf plot a good first step in guessing the order q of the process.

• An \underline{autoregressive} process of order p, denoted AR(p) has the form where again z_{t} is a mean-zero white noise process. In terms of the backshift operator, By_{t} = y_{t−1}, we may write

(1 − α_{1}B − α_{2}B^{2} − ⋯ − α_{p}B^{p})y_{t} = z_{t}.

If E[y_{t}] = μ ≠ 0, then we would write

(1 − α_{1}B − α_{2}B^{2} − ⋯ − α_{p}B^{p})(y_{t} − μ) = z_{t}.

Not all values of the coefficients α_{k} lead to a stationary series. For example we saw last time that the random walk y_{t} = y_{t−1} +z_{t} was not a stationary process. More generally the complex roots of the polynomial

Φ(x) = 1 − α_{1}x − α_{2}x^{2} − ⋯ − α_{p}x^{p}

must have all p of its roots outside the unit circle in the complex plane. In the case of an AR(1) process, this is easy to show!

One can show that the autocovariance function for a stationary AR(p) process satisfies the “Yule-Walker” relations

∑_{k=1}^{p} γ(k − i)α_{k} = γ(k) i = 1, 2, …, p,

We can use the Yule-Walker equations to solve for γ(k) in terms of the α’s, if they are known, or solve for the α’s in terms of the γ(k)’s. This suggests a simple way to estimate the coefficients of an AR process. Also,

σ_{y}^{2} = ∑_{k=1}^{p} γ(k)α_{k} + σ_{z}^{2}.

As an exercise, show that the autocorrelation function of an AR(1) process y_{t} = αy_{t−1} + z_{t} satisfies ρ_{k} = αρ_{k−1}, and so, inserting ρ_{0} = 1, we get ρ_{k} = α^{k}. What does ρ_{k} look like for an AR(2) process? Here are plots of simulated AR(1) and AR(2) processes.

The acf plot of an AR(1) decays exponentially to zero, possibly with some oscillation in sign. The AR(2) acf plot also decays to zero (this gradual decay in the acf plot distinguishes pure-AR from pure-MA processes) but it isn’t easy to tell from the plot what the order p of the process should be.

Partial Autocorrelation Function

Here are some more examples of AR processes and their acf’s:

• The following graph shows plots of two AR processes,

y_{t} = .7y_{t−1} + .2y_{t−2} + z_{t}

and

y_{t} = .7y_{t−1} + z_{t}

and their autocorrelation functions. Notice that the ACF’s have long lags and it is not easy to diagnose anything.

402 > tsplot(dd <- arima.sim(200, model=list(ar=c(.7, .2))))
402 > acf(dd)
402 > tsplot(dd <- arima.sim(200, model=list(ar=c(.7))))
402 > acf(dd)

• Next we plot two more MA processes and their acf’s:

402 > tsplot(dd <- arima.sim(200, model=list(ma=c(.6))))
402 > acf(dd)
402 > tsplot(dd <- arima.sim(200, model=list(ma=c(.6, .4))))
402 > acf(dd)

The two models that we fit are

y_{t} = .6z_{t−1} + z_{t}

and

y_{t} = .6z_{t−1} + .4z_{t−2} + z_{t}

Just as before, that the ACF tails off after lag 1 and 2. That is the ACF provides a diagnostic of the order of the MA process.

A simple (too simple!) conclusion is that if the ACF is zero beyond a certain lag, then that is evidence for an MA process of the appropriate order. That turns out to be a useful observation, but it doesn’t help us with AR processes.

The form of an autoregressive model, suggests that instead of plotting the covariances between lagged observations, we should be plotting the coefficients α_{k} themselves: essentially, we want to fit AR models of orders 1, 2, 3, and so on. Then at lag 1 we plot the coefficient α_{1} of the fitted AR(1) model. At lag 2 we plot the coefficient α_{2} of the AR(2) model, and so on. When p is too high, the plotted coefficients should be essentially zero. This type of plot is called a partial autocorrelation function (PACF), and it works really well.

To actually calculate the PACF (and plot it in SPLUS) requires the use of the acf function with an extra argument.

acf(dd,type="p")

I have produced four plots, two of AR(2) processes and two of AR(1) processes. Notice that these plots look like the ACF plots of MA processes.

402 > acf(arima.sim(200, model=list(ar=c(.7, .2))), type="p")
402 > acf(arima.sim(200, model=list(ar=c(.7, .2))), type="p")
402 > acf(arima.sim(200, model=list(ar=c(.7))), type="p")
402 > acf(arima.sim(200, model=list(ar=c(.7))), type="p")

Maybe the PACF of a MA process will look like the ACF of an AR process? This sort of duality would be quite useful. Lets look and see. The simulations were generated by

402 > acf( arima.sim(200, model=list(ma=c(-.9))), type="p")

402 > acf( arima.sim(200, model=list(ma=c(.9))), type="p")
402 > acf( arima.sim(200, model=list(ma=c(.9, .6))), type="p")
402 > acf( arima.sim(200, model=list(ma=c(.9, .6))), type="p")

The following plots really do look like the ACF of AR’s. There is a mathematical duality between ACF’s and PACF’s, but the important thing to note is that each of them is quite good at helping to diagnose the order of a pure MA or a pure AR process.

Thus,

• To get an idea what the order q of a pure MA(q) process should be, use an ACF plot.

• To get an idea what the order p of a pure AR(p) process should be, use a PACF plot.

There are better ways to choose p and q, but these make useful illustrations and can help you diagnose whether you are “way off”.