\[ \newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]} \newcommand{\Var}[1]{\mathrm{Var}\left[ #1 \right]} \newcommand{\SampleVar}[1]{\widehat{\mathrm{Var}}\left[ #1 \right]} \newcommand{\Cov}[1]{\mathrm{Cov}\left[ #1 \right]} \newcommand{\TrueRegFunc}{\mu} \newcommand{\EstRegFunc}{\widehat{\TrueRegFunc}} \DeclareMathOperator{\tr}{tr} \DeclareMathOperator*{\argmin}{argmin} \DeclareMathOperator{\det}{det} \newcommand{\TrueNoise}{\epsilon} \newcommand{\EstNoise}{\widehat{\TrueNoise}} \newcommand{\SignalNoise}{N} \]

In our previous episodes: Decomposing data into components

• Eigenvectors of the smoother
  • what’s kept by smoothing?
  • what shows up in the residuals?
• Eigenvectors of the variance matrix for multivariate data
  • what are the principal components?
  • == what are the directions of maximum variance?
  • == what are the uncorrelated components?
  • == what components should we keep when reducing dimension?

Today: Decomposing data into periodic components

• Helps with periodic trends/seasonality
• Also helps with covariance estimation
• And with any sort of linear operation

Sine waves

We get data at times \(1, 2, \ldots n\)

curve(sin(2 * pi * x/n), xlim = c(0, n), ylim = c(-2, 2), type = "p", pch = 16, 
    cex = 0.5, col = "blue")

Sine waves

curve(sin(2 * pi * x/n), xlim = c(0, n), ylim = c(-2, 2), type = "p", pch = 16, 
    cex = 0.5, col = "blue")
curve(sin(2 * pi * 3 * x/n), type = "b", pch = 16, cex = 0.5, col = "yellow3", 
    add = TRUE)

Sine waves

curve(sin(2 * pi * x/n), xlim = c(0, n), ylim = c(-2, 2), type = "p", pch = 16, 
    cex = 0.5, col = "blue")
curve(sin(2 * pi * 3 * x/n), type = "p", pch = 16, cex = 0.5, col = "yellow3", 
    add = TRUE)
curve(sin(2 * pi * x/n) + sin(2 * pi * 3 * x/n), type = "p", pch = 16, cex = 0.5, 
    col = "green", add = TRUE)

An ambition

• Examples like this suggest a new project: decompose arbitrary time series into a sum of sine waves
  • i.e., write anything, periodic or not, as a sum of periodic components

Sines average out to zero

For any integer \(\omega\), \[\begin{eqnarray} \sum_{t=0}^{n-1}{\sin{2\pi\omega \frac{t}{n}}} & = & 0\\ \end{eqnarray}\] because \[ \sin{2\pi\omega \frac{t}{n}} = -\sin{2\pi\omega \frac{t+n/2\omega}{n}} \] i.e., time points come in pairs which cancel each other out in the sum

• ditto for averages of cosines \[ \sum_{t=0}^{n-1}{\cos{2\pi\omega\frac{t}{n}}} = 0 ~ \text{unless} ~ \omega =0 \]

Sines of different frequencies are orthogonal

• The inner product of two sine waves is usually 0: \[\begin{eqnarray} \sum_{t=1}^{n}{\sin{\left(2\pi\omega_1 \frac{t}{n}\right)} \sin{\left(2\pi\omega_2 \frac{t}{n}\right)}} & = & \frac{1}{2}\sum_{t=1}^{n}{\cos{\left(2\pi(\omega_1-\omega_2) \frac{t}{n}\right)} - \cos{\left(2\pi(\omega_1+\omega_2)\frac{t}{n}\right)}}\\ & = & 0 ~ \text{unless} ~ \omega_1=\omega_2 ~ \text{or}~ \omega_1 = -\omega_2 \end{eqnarray}\]

  • If \(\omega_1 = \omega_2 \neq 0\), we get \(n/2\)
  • If \(\omega_1 = -\omega_2 \neq 0\), we get \(-n/2\)
  • If \(\omega_1 = \omega_2 = 0\), we get 0

• ditto for inner products of cosines

We can’t resolve either very high or very low frequency sine waves

• High frequency: If \(\sin{\phi + 2\pi\omega t/n}\) goes through a complete cycle between \(t\) and \(t+1\), then it looks exactly like \(\sin{\phi}\) everywhere
• Low frequency: If \(\omega\) is too small, we don’t see even half a cycle, and we can’t tell the difference between that and a constant plus higher-frequency terms

Orthogonal vectors give us a basis

• Remember from linear algebra: vectors \(\vec{v}_1, \vec{v}_2, \ldots \vec{v}_n\) are a basis when, for any vector \(\vec{x}\), \[ \vec{x} = \sum_{i=1}^{n}{c_i \vec{v}_i} \]
• In an \(n\)-dimensional space, any set of \(n\) orthogonal vectors is a basis
  • Also true of any \(n\) linearly independent vectors
• So we can use the sines and cosines as a basis
  • Not normalized, but that’s OK

Re-writing the data in terms of sines and cosines

• Re-write any vector in terms of the basis: \[ x(t) = \sum_{\omega=0}^{n/2-1}{a_{\omega}\sin{\left(2\pi\omega \frac{t}{n}\right)} + b_{\omega}\cos{\left(2\pi\omega\frac{t}{n}\right)}} \]

• Why only up to \(n/2-1\)?
• \(b_0 = \overline{x}\) (why?) and \(a_0=0\) (why?)

• The \(a\) and \(b\) coefficients are linked…

Re-writing the data in terms of sines

Use the trig identities to eliminate cosines in favor of sines:

\[ x(t) = \sum_{\omega=0}^{n/2-1}{\alpha_{\omega} \sin{\left(2\pi\frac{t}{n} + \phi_{\omega}\right)}} \] with \[ \alpha_{\omega} = \sqrt{a_{\omega}^2+b_{\omega}^2} \]

Re-writing the data in terms of complex exponentials

• \(e^{i \theta} = \cos{\theta} + i \sin{\theta}\)
• Re-write to eliminate \(\sin\) in favor of exponentials: \[ x(t) = \sum_{\omega=0}^{n-1}{c_{\omega} e^{2\pi i \omega t/n}} \]

The Fourier transform

• The Fourier transform of \(x(t)\), \(t=0, \ldots n-1\), is \[ \tilde{x}(\omega) \equiv \sum_{t=0}^{n-1}{x(t) e^{-2\pi i t \omega/n}} \]
• These numbers are the Fourier coefficients
• \(\tilde{x}(\omega)= r e^{i\phi}\) with \(r, \phi\) real
  • \(r =\) amplitude
  • \(\phi =\) phase
• This is the frequency domain

The inverse Fourier transform

• The inverse Fourier transform recovers \(x(t)\) from \(\tilde{x}(\omega)\): \[ x(t) = \frac{1}{n}\sum_{\omega=0}^{n/2-1}{\tilde{x}(\omega) e^{2\pi i t \omega/n}} \]

• We are back in the time domain
  • Always comes out real
• The un-normalized inverse transform is \[ \sum_{\omega=0}^{n/2-1}{\tilde{x}(\omega) e^{2\pi i t \omega/n}} \]

An example

n <- 60
x <- sin(2 * pi * (1:60)/n) + 0.5 * sin(2 * pi * 3 * (1:60)/n)
plot(x, xlab = "t", ylab = expression(x(t)))

An example

An example

In R

• You could do everything directly from the definitions, but this is slow
• Use the fast Fourier transform (FFT), from the fft function
  • The trick is cute, but it’s just a computing trick
  • The inverse=TRUE option gives the unnormalized inverse
• Because R starts counting from 1, not 0, give vector x of length n,

fft(x)[h] = sum_{k=1}^n {x[k]*exp(-2*pi*1i*(k-1)*(h-1)/n)}

so entry h in the output of fft contains the Fourier coefficient for frequency (h-1)/n

In R

plot(x, pch = 16)

In R

plot(x, pch = 16)
points(Re(fft(fft(x), inverse = TRUE)/length(x)), col = "red", pch = 25)

Why do we care?

• Periodic components show up!
• Differentiation and integration are easy: \[ \frac{d}{dt}\left( c_{\omega} e^{2\pi i \omega t/n}\right) = \frac{2\pi i \omega}{n} c_{\omega} e^{2\pi i \omega t/n} \] so \[ \widetilde{x^{\prime}}(\omega) = \frac{2\pi i \omega}{n} \tilde{x}(\omega) \]
  • (similarly for integration)
• Shifts in time are easy: if \(y(t) = x(t-t_0)\) then \[ \tilde{y}(\omega) = \tilde{x}(\omega) e^{-2\pi i t_0 \omega/n} \]
• Convolutions are easy

Convolutions

• Fix any function \(h\) of time
• The convolution of \(x(t)\) with \(h\) is a new function of time: \[\begin{eqnarray} (x \star h)(t) &\equiv & \sum_{s=0}^{n-1}{x(s) h(t-s)} \end{eqnarray}\]
  • A weighted sum of \(x\) values around \(t\), with weights given by \(h\)
  • We’ve seen this sort of thing before!
  • Moving averages = convolution with indicator function
  • Kernel smoothing = convolution with the kernel

Convolutions and Fourier transforms

• Convolution in the time domain is multiplication in the frequency domain: \[ \widetilde{x \star h}(\omega) = \tilde{x}(\omega) \tilde{h}(\omega) \]
  • Reverse is also true: convolution in the frequency domain is multiplication in the time domain
• Doing convolution quickly:
  • Take the Fourier transforms of \(x\) and \(h\) (using fft)
  • Multiply the Fourier coefficients
  • Take the inverse transform
  • Implemented in the convolve function

Autocovariance and Fourier transforms

• Suppose \(X(t)\) stationary, with autocovariance \(\gamma(h)\)
• We can do a Fourier transform on \(\gamma\): \[ \tilde{\gamma}(\omega) = \sum_{h=0}^{n-1}{\gamma(h) e^{-2\pi i h \omega/n}} \]
• How does this relate to \(\tilde{x}\)?

Power spectrum

• \(X(t)\) is random, so \(\tilde{X}(\omega)\) is also random
• Amplitude \(|\tilde{X}(\omega)|\) is random
• Power \(|\tilde{X}(\omega)|^2\) is random
• Power spectrum \[ S(\omega) = \Expect{|\tilde{X}(\omega)|^2} \] is not random

The Wiener-Khinchin theorem:

• The power spectrum is the Fourier transform of the autocovariance: \[ \tilde{\gamma}(\omega) = S(\omega) = \Expect{|\tilde{X}(\omega)|^2} \]
• Uses:
  • If we can estimate one side, we can estimate the other
  • If we can estimate \(\Expect{|\tilde{X}(\omega)|^2}\), we can estimate \(\gamma(h)\)

The periodogram

• The plot of \(|\tilde{x}(\omega)|^2\) against \(\omega\) is the periodogram
  • As \(n \rightarrow \infty\), \(\Expect{|\tilde{x}(\omega)|^2} \rightarrow S(\omega)\)
  • but \(\Var{|\tilde{x}(\omega)|^2} \not\rightarrow 0\)
• Solution: smooth the periodogram
• Basic R command: spectrum
  • Some sensible defaults about smoothing
  • annoyingly, different normalization of frequencies than fft

In R

spectrum(x)

In R

library(gstat)
data(wind)
wind <- wind[!(wind$month == 2 & wind$day == 29), ]
spectrum(wind[, "DUB"], span = 10)

In R

library(gstat)
data(wind)
wind <- wind[!(wind$month == 2 & wind$day == 29), ]
spectrum(wind[, "DUB"], span = 10)
abline(v = 1/365)

Summary

Some cautions:

• Lots of variants in notation/definitions
  • some people absorb the \(2\pi\) into the definition of \(\omega\)
  • some people absorb the \(1/n\) into the definition of \(\omega\)
  • some people write \(k\) or \(h\) or \(f\) instead of \(\omega\)
  • some people normalize the basis functions

Backup: Eliminating sine+cosine in favor of sine

Backup: why \(x(t) e^{-2\pi i \omega t/n}\)?

The hand-wavy explanation:

• When we take the inner product of two real vectors \(\vec{u}\) and \(\vec{v}\), it’s the same as \(\vec{u}^T \vec{v}\) (\([1\times p]\) \([p\times 1]\))
• But the inner product of anything with itself is a (squared) length, so \(\geq 0\)
• If \(\vec{v} = a+bi\), then \(\vec{v}^T \vec{v} = (a+bi)(a+bi) = a^2 - b^2 + 2iab\), which doesn’t make sense as a length
• Solution: don’t just transpose one vector, transpose it and take its complex conjugate
  • Sometimes called finding the adjoint vector
  • Often written \(\vec{v}^*\) or \(\vec{v}^\dagger\)
• Now \(\vec{v}^* \vec{v} = (a-bi)(a+bi) = a^2 + b^2\), which makes sense as a length
• We want to take the inner product of the vector of \(x(t)\)’s with the basis vector we get from \(e^{2\pi i \omega t/n}\), so we need to complex-conjugate the basis vector

Backup: why \(x(t) e^{-2\pi i \omega t/n}\)?

The “shut up and calculate” explanation:

• The (complex) numbers \(\tilde{x}(\omega)\) are defined to be the ones which give us back the (real) numbers \(x(t)\): \[ x(t) = \frac{1}{n}\sum_{\omega=0}^{n/2-1}{\tilde{x}(\omega) e^{2\pi i t \omega/n}} \]
• Pick your favorite frequency \(\nu\), and multiply both sides by \(e^{-2\pi i t \nu/n}\): \[ x(t) e^{-2\pi i t \nu/n} = \frac{1}{n}\sum_{\omega=0}^{n/2-1}{\tilde{x}(\omega) e^{2\pi i t (\omega-\nu)/n}} \]
• Sum over all \(t\): \[ \sum_{t=1}^{n}{x(t) e^{-2\pi i t \nu/n}} = \frac{1}{n}\sum_{t=1}^{n}{\sum_{\omega=0}^{n/2-1}{\tilde{x}(\omega) e^{2\pi i t (\omega-\nu)/n}}} \]
• Exchange the order of summation on the right-hand side: \[ \sum_{t=1}^{n}{x(t) e^{-2\pi i t \nu/n}} = \frac{1}{n}\sum_{\omega=0}^{n/2-1}{\tilde{x}(\omega)\sum_{t=1}^{n}{ e^{2\pi i t (\omega-\nu)/n}}} \]
• Now remember that sines and cosines average to zero, except at zero frequency: \[ \sum_{t=1}^{n}{x(t) e^{-2\pi i t \nu/n}} = \frac{1}{n}\tilde{x}(\nu)n = \tilde{x}(\nu) \]