Start with the Karhunen-Loeve transform

• Given: \(X(t)\) defined on a finite domain (e.g., interval) \(D\), \(\Expect{X(t)}=0\) (centering)
  • Substitute \(\mathbf{r}\) for \(t\), or \((\mathbf{r}, t)\), etc.
• Desired \(X(t) = \sum_{i=1}^{\infty}{S_i \phi_i(t)}\) with random coefficients \(S_i\)
  • Implies \(\Expect{S_i}=0\)
• Require orthonormality of the functions: \(\int_{D}{\phi_i(t) \phi_j(t) dt} = \delta_{ij}\)
• Implies \(\Cov{X(t), X(s)} = \sum_{i=1}^{\infty}{\sum_{j=1}^{\infty}}{\phi_i(t) \phi_j(s) \Cov{S_i, S_j}}\)
• Require minimal mean squared error, so we pick the \(\phi_i\) to minimize \(\Expect{\int_{D}{\left(X(t) - \sum_{i=1}^{q}{S_i \phi_i(t)}\right)^2 dt}}\) at each \(q\), subject to the orthonormality constraint
• Implies \(\int_{D}{\Cov{X(t), X(s)} \phi_i(s) ds} = \lambda_i \phi_i(t)\)
  • An eigenproblem again, but now solved by an eigenfunction, not an eigenvector
  • Eigenproblems for integral operators are classic early-20th-century math (Courant and Hilbert 1953, ch. III)
    • \(\Cov{X(t), X(s)}\) is the kernel of the integral operator (and yes, stuff like this is where our “kernel” comes from)
• Implies \(\Cov{S_i, S_j} = \lambda_i \delta_{ij}\)
  • \(\therefore\) \(\Cov{X(t), X(s)} = \sum_{i}^{\infty}{\phi_i(t) \phi_i(s) \lambda_i}\)
• We need some sort of estimate of the covariance function…

Some things to notice about Karhunen-Loeve / PCA

• We are transforming to new coordinates
• The expected value of the transform is zero
  • Because we centered \(X\); otherwise, add the constant function to soak up the mean
• The variance of the transformed coordinates is not zero
  • We’re dealing with a random process
• The variance of transformed coordinates tells us a lot about the underlying process
  • \(\Var{S_i} = \lambda_i\) and \(\lambda_i\) is an eigenvalue
• The covariance across the transformed coordinates is zero
• The transformation is different for every process

Let’s use sine waves!

• Let’s fix the basis functions, independent of the random process
• There are many, many choices
• Sine waves are particularly use for stationary processes
• Stationary = the distribution does not change over time
• Weakly stationary = second-order stationary = mean, variance and covariance do not change over time
• This plays nicely with sine waves because \(e^{2\pi i \omega (t-t_0)} = e^{-2\pi i \omega t_0} e^{2\pi i \omega t}\)
  • Shifting in time changes the phase of the wave, in a trivial way, but nothing else
  • Most fixed choices of basis function are not so nice!
• “Decompose an arbitrary function into a linear combination of sine waves” = Fourier analysis
• Pros: the math is nice, and well-developed
• Cons:
  • We’ll need a lot of coefficients / transformed coordinates
  • The new coefficients won’t be uncorrelated (in general)
  • We need complex-valued coefficients

The Fourier transform (in 1D)

• For continuous time, \[ \tilde{x}(\omega) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}{x(t) e^{2\pi i \omega t} dt} \]
  • Note: a continuous infinity of basis functions
• For discrete time, \(j=1, \ldots n\), \[ \tilde{x}_k = \frac{1}{\sqrt{n}}\sum_{j=1}^{n}{x_j e^{2\pi i k j / n}} \]
  • Only \(n\) (orthonormal) basis functions
  • If the physical time between measurements (\(\equiv\) sampling interval) is \(\Delta t\), then the corresponding physical frequencies \(\omega_k = \frac{k}{n \Delta t}\)
• The Fourier transform is a linear operator say \(\mathcal{F}\)
  • Maps functions on the time domain into functions on the frequency domain
    • Strictly speaking \(\mathcal{F}\) maps complex-valued functions of the real line to complex-valued functions of the real line, \(\mathcal{F}: (\mathbf{C}^{\mathbf{R}}) \mapsto (\mathbf{C}^{\mathbf{R}})\), but domain and the range mean different things

The Fourier transform (in 1D)

• The inverse Fourier transform is \[ \tilde{\tilde{x}}(t) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}{\tilde{x}(\omega) e^{-2\pi i \omega t} d\omega} = x(t) \] or \[ \tilde{\tilde{x}}_j = \frac{1}{\sqrt{n}}\sum_{k=-n/2}^{n/2}{\tilde{x}_k e^{-2\pi i k j/n}} = x_j \]
  • This is, again, a linear operator
  • The inverse transform takes frequency-domain functions to time-domain functions
  • The inverse transform looks almost just like the forward transform
  • WARNING: Some people use slightly different conventions for forward/inverse FT; some people absorb the \(2\pi\) into \(\omega\) or \(k\)

Convolution

• The convolution of \(f\) and \(g\) is a new function, \(f*g\): \[ (f*g)(t) \equiv \int{f(s) g(t-s) ds} = (g*f)(t) \]
  • EXERCISE: Prove \(f*g = g*f\)
  • similarly in discrete time
• If \(g\) is a pdf, then \((f*g)\) is a weighted average of \(f\) at each point
• Convolving \(f\) with a pdf gives a new, smoother function
  • Kernel (= Nadaraya-Watson) smoothing is convolving with a pdf (= the kernel)
• If \(\mathrm{pdf}(X) = f\) and \(\mathrm{pdf}(Z) = g\) and \(X\) and \(Z\) are independent, then \(Y=X+Z\) has pdf \(f*g\)
• The Fourier transform turns convolution into multiplication: \[ \widetilde{f*g}(\omega) = \tilde{f}(\omega) \tilde{g}(\omega) \]
  • Conversely, it turns multiplication into convolution: \(\widetilde{fg} = \tilde{f} * \tilde{g}\)
  • Proving this is left as an EXERCISE
    • See the reading if you get stuck, but try it on your own

Estimating the power spectrum

The periodogram

• Take the discrete Fourier transform of \(X_j\), get \(\tilde{X}_k\)
• Set \(\hat{S}(\omega_k) = |\tilde{X}_{k}|^2\)
• This is unbiased
• And \(\Cov{\hat{S}(\omega), \hat{S}(\nu)} \rightarrow 0\) as \(n\rightarrow \infty\)
• But we have limited frequency resolution \(\Rightarrow\) “leakage” of power from nearby frequencies
• And \(\Var{\hat{S}(\omega)} \not\rightarrow 0\) as \(n\rightarrow \infty\)

Practically

• Smooth the periodogram by convolving it with a pdf
  • Deliberately add bias, to reduce variance
• This is actually where statisticians first used kernel methods, in the 1950s
  • The people who did this went on to invent kernel density estimation (Rosenblatt 1956; Parzen 1962)
• Smoothing less and less as we get more data gives consistent estimates of the power spectrum

Filtering, smoothing, trends

• The inverse Fourier transform recovers the original time series exactly
• We can deliberately change some of the Fourier coefficients, and then the inverse transform will give us something else
• One common trick: \(\tilde{x}(\omega) \rightarrow \tilde{x}(\omega) \kappa(\omega)\) where \(\kappa(u) \geq 0\), \(\kappa(u) \rightarrow 0\) as \(|u|\rightarrow \infty\)
  • \(\kappa\) is sometimes, but not always, a pdf
• Then \(\hat{x} = \mathcal{F}^{-1} (\tilde{x} \kappa)\)
• The fact that \(\kappa(u) \rightarrow 0\) as \(|u| \rightarrow \infty\) makes this a low-pass filter
  • Low-frequency components are preserved, high-frequency components are attenuated
  • What would a high-pass filter be?
• Multiplication in the frequency domain is convolution in the time domain: \(\hat{x} = x * \tilde{\kappa}\)
  • This is a smoother again
• One way to estimate a trend

Non-Stationarity

1. We can de-trend, and then apply Fourier analysis to the fluctuations.
2. We can break the data up into segments, and hope it’s (approximately) stationary within each segment.
   1. A spectrogram is a plot showing (local) Fourier transforms for successive time segments
   2. This leads to a test of stationarity: we should get the same Fourier spectrum from each segment, plus or minus sampling error.
      1. Fuentes, Guttorp, and Sampson (2007), sec. 1.6.2 gives one version of this idea, with conditions under which we can make it into an ANOVA test and use the standard machinery
      2. Kantz and Schreiber (2004), sec. 2.2–2.3 and 13.2, give what is (essentially) a bootstrap version of the same idea
3. Give up on the Fourier basis in favor of wavelets
   • Fourier basis functions have equal amplitude everywhere
   • Come up with wiggly, wave-like functions that are localized
   • It’s just another set of basis functions, so we just take inner products again

Backup: Simulating from a power spectrum

Theorem: Suppose \(g(\omega) \geq 0\), and define \(G(\omega) = \int_{-\infty}^{\omega}{g(\nu) d\nu}\). Further suppose that \(G(-\infty) = 0\) and \(G(\infty)=\sigma^2/\sqrt{2\pi} > 0\). Then there exists a centered, weakly-stationary stochastic process \(X\) with autocovariance function \(\tilde{g}\). In particular, \(\Var{X(t)} = \sigma^2\) for all \(t\), and \(\Cov{X(t_1), X(t_2)} = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}{ e^{-2\pi i (t_2-t_1) \omega} g(\omega) d\omega}\).

Sketch proof: Draw \(W \sim G/(\sigma^2 \sqrt{2\pi})\). That is, the pdf of \(W\) is \(g(\omega)/(\sigma^2 \sqrt{2 \pi})\), and \(\Pr{\left(W leq \omega\right)} = G(\omega)/(\sigma^2 \sqrt{2\pi})\). Draw \(A \sim \mathrm{Unif}(0, 2\pi)\), independent of \(N\). Set \(X(t) \equiv \frac{\sigma} e^{iA} e^{-iNt}\). It’s easy to check that \(\Expect{e^{iA}} = 0\), so \(\Expect{X(t)} = 0\) too. And \[ \Expect{\overline{X(t_1)} X(t_2)} = \sigma^2 \Expect{e^{iN(t_2 -t_1)}} = \int_{-\infty}^{\infty}{e^{-2\pi i \omega (t_2-t_1)} \frac{g(\omega)}{\sqrt{2\pi}} d\omega} \] as desired.

Notes:

• This is one of those places where things would be cleaner if I didn’t have the factor of \(1/\sqrt{2\pi}\) in the inverse Fourier transform
• We need to complex-conjugate either \(X(t_1)\) or \(X(t_2)\) to get a covariance function for complex-valued processes.
• You can modify this to handle dis-continuous power spectra, too.

Backup: Simulating from a power spectrum (cont’d.)

• The construction on the previous slide might not be very satisfying, because every realization is just a sine wave, and the stochasticity only becomes apparent across realizations
• Here’s a related idea: for each \(\omega\), \(A(\omega) \mathrm{Unif}(0, 2\pi)\), IIDly, except that \(A(-\omega) = \overline{A(\omega)}\) \[ X(t) = (2\pi)^{-1/4} \int{e^{iA(\omega)} e^{-2\pi i \omega t} \sqrt{g(\omega)} d\omega} \] Now \[ \Expect{\overline{X(t_1)} X(t_2)} = \frac{1}{\sqrt{2\pi}} \int{\left( \int{\Expect{e^{-iA(\omega)} e^{iA(\nu)}} e^{-2\pi i (\omega t_i - \nu t_2)} \sqrt{g(\nu)} d\nu} \right) \sqrt{g(\omega)} d\omega} = \frac{1}{\sqrt{2\pi}} \int{e^{-2\pi i \omega (t_2 - t_1)} g(\omega) d\omega} \] (How did I get rid of the double integral?)