\[ \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{\Signal}{S} \newcommand{\SignalNoise}{N} \newcommand{\AutoCov}{\gamma} \newcommand{\NoiseAutoCov}{\xi} \newcommand{\LinearExpect}[1]{\mathbb{L}\left[ #1 \right]} \newcommand{\Periodicity}{\tau} \]

In our previous episodes

General approach to optimal linear prediction
- Predict \(Y\) from \(\vec{Z} = [Z_1, Z_2, \ldots Z_p]\)
- Best linear predictor \(\LinearExpect{Y|\vec{Z}} = \alpha + \vec{\beta} \cdot \vec{Z}\)
- Best \(\alpha = \Expect{Y} - \vec{\beta} \cdot \Expect{\vec{Z}}\)
- Best \(\vec{\beta} = \Var{\vec{Z}}^{-1} \Cov{\vec{Z}, Y}\)
Previously:
- Interpolating and extrapolating time series
- Interpolating and extrapolating random fields
Today:
- Extracting a random signal from random noise
- Extracting a deterministic, periodic trend from noise

Optimal linear filtering for time series

Given: \(X(t_1), X(t_2), \ldots X(t_n)\), for short \(X_1, \ldots X_n\)
Assumption: \(X(t) = \Signal(t) + \SignalNoise(t)\), with \(\Expect{\SignalNoise(t)} = 0\)
- Signal \(\Signal(t)\) plus noise \(\SignalNoise(t)\)
- NOTE: We’re assuming signal \(S\) is also a random process!
Desired: prediction of \(\Signal(t_0)\), for short \(\Signal_0\)
- This is filtering
- \(t_0\) could be one of the \(t_i\), or between them, or into the future or past

Optimal linear filtering for time series

\[\begin{eqnarray} \hat{\Signal}(t_0) & = & \alpha + \vec{\beta} \cdot \left[\begin{array}{c} X_1 \\ X_2 \\ \vdots \\ X_n \end{array}\right]\\ \alpha & = & \Expect{\Signal_0} - \vec{\beta} \cdot \left[\begin{array}{c} \Expect{X_1}\\ \Expect{X_2} \\ \vdots \\ \Expect{X_n}\end{array}\right]\\ \vec{\beta} & = & {\left[\begin{array}{cccc} \Var{X_1} & \Cov{X_1, X_2} & \ldots & \Cov{X_1, X_n}\\ \Cov{X_1, X_2} & \Var{X_2} & \ldots & \Cov{X_2, X_n}\\ \vdots & \vdots & \ldots & \vdots\\ \Cov{X_1, X_n} & \Cov{X_2, X_n} & \ldots & \Var{X_n}\end{array}\right]}^{-1} \left[\begin{array}{c} \Cov{\Signal_0, X_1}\\ \Cov{\Signal_0, X_2}\\ \vdots \\ \Cov{\Signal_0, X_n}\end{array}\right] \end{eqnarray}\]

This should look familiar

It’s almost the same as linear prediction for time series
It’s exactly the same as linear prediction of one time series from another

Simple case, to build intuition

Assume \(\Signal\) is stationary, with \(\Expect{S(t)} = 0\), autocovariance function \(\AutoCov(h)\)
Assume \(\SignalNoise\) is stationary, with variance \(\sigma^2\) and no autocorrelation or correlation with \(S\)

\[\begin{eqnarray} \Cov{X(t), X(t+h)} & = & \Cov{\Signal(t)+\SignalNoise(t), \Signal(t+h)+\SignalNoise(t+h)}\\ & = & \Cov{\Signal(t), \Signal(t+h)} + \Cov{\SignalNoise(t), \SignalNoise(t+h)}\\ & & + \Cov{\Signal(t), \SignalNoise(t+h)} + \Cov{\Signal(t+h), \SignalNoise(t)}\\ & = & \AutoCov(h) + \sigma^2\delta_{h0}\\ \Cov{\Signal(t), X(t+h)} & = & \Cov{\Signal(t), \Signal(t+h) + \SignalNoise(t+h)}\\ & = & \Cov{\Signal(t), \Signal(t+h)} + \Cov{\Signal(t), \SignalNoise(t+h)}\\ & = & \AutoCov(h) \end{eqnarray}\]

(Last time I said it’s common to add an extra “nugget” of auto-covariance at lag 0 for the observations; that’s handling measurement noise)

One observation, one estimate

Given: \(X(t)\)
Desired: estimate/prediction of \(\Signal(t)\)

\[\begin{eqnarray} \hat{\Signal}(t) & = & \beta X(t)\\ \beta & = & \frac{\Cov{\Signal(t), X(t)}}{\Var{X(t)}}\\ &= & \frac{\AutoCov(0)}{\AutoCov(0) + \sigma^2} < 1\\ \hat{\Signal}(t) & = & \frac{\AutoCov(0)}{\AutoCov(0) + \sigma^2} X(t) \end{eqnarray}\]

\[ \frac{\AutoCov(0)}{\AutoCov(0) + \sigma^2} = \frac{1}{1+\sigma^2/\AutoCov(0)} \] \(\therefore\) everything turns on the ratio of noise variance (\(\sigma^2\)) to signal variance (\(\AutoCov(0)\)), a.k.a. signal-to-noise ratio

Two observations, one estimate

Given: \(X(t-1)\), \(X(t)\)
Desired: estimate/prediction of \(\Signal(t)\)
Estimate is \[\begin{eqnarray} \hat{\Signal}(t) & = & \beta_0 X(t) + \beta_1 X(t-1)\\ \left[ \begin{array}{cc} \beta_0 \\ \beta_1 \end{array}\right] & = & \left[\begin{array}{cc} \Var{X(t)} & \Cov{X(t), X(t-1)}\\ \Cov{X(t-1), X(t)} & \Var{X(t-1)} \end{array}\right]^{-1} \left[\begin{array}{c} \Cov{X(t), \Signal(t)} \\ \Cov{X(t-1), \Signal(t)}\end{array} \right] \end{eqnarray}\]

Two observations, one estimate

\[\begin{eqnarray} \left[ \begin{array}{cc} \beta_0 \\ \beta_1 \end{array}\right] & = & \left[\begin{array}{cc} \Var{X(t)} & \Cov{X(t), X(t-1)}\\ \Cov{X(t-1), X(t)} & \Var{X(t-1)}\end{array}\right]^{-1} \left[\begin{array}{c} \Cov{X(t), \Signal(t)} \\ \Cov{X(t-1), \Signal(t)}\end{array} \right]\\ & = & \left[\begin{array}{cc} \AutoCov(0)+\sigma^2 & \AutoCov(1)\\ \AutoCov(1) & \AutoCov(0)+\sigma^2\end{array}\right]^{-1} \left[\begin{array}{c} \AutoCov(0) \\ \AutoCov(1) \end{array} \right]\\ &=& \frac{1}{(\AutoCov(0)+\sigma^2)^2-\AutoCov^2(1)}\left[\begin{array}{cc} \AutoCov(0)+\sigma^2 & -\AutoCov(1)\\ -\AutoCov(1) & \AutoCov(0)+\sigma^2\end{array}\right]\left[\begin{array}{c} \AutoCov(0) \\ \AutoCov(1) \end{array} \right]\\ & = & \frac{1}{(\AutoCov(0)+\sigma^2)^2-\AutoCov^2(1)}\left[\begin{array}{c} (\AutoCov(0)+\sigma^2)\AutoCov(0) - \AutoCov^2(1) \\ \AutoCov(1)\sigma^2\end{array}\right] \end{eqnarray}\]

What’s going on here?

Slope on \(X(t) < 1\) \(\Rightarrow\) Estimate of \(S(t)\) should not change 1-for-1 with \(X(t)\)
Weak noise \(\Rightarrow\) pay more attention to the data
Strong correlation in the signal \(\Rightarrow\) pay more attention to \(X(t-1)\) when estimating \(\Signal(t)\)
- Reduces to the one observation, one signal case when \(\AutoCov(1) = 0\)

The general pattern: the Wiener filter

Assume \(X(t) = \Signal(t) + \SignalNoise(t)\)
Assume \(\SignalNoise\) uncorrelated with \(\Signal\)
Implications: \[\begin{eqnarray} \Cov{X(t), X(t+h)} & = & \Cov{\Signal(t), \Signal(t+h)} + \Cov{\SignalNoise(t), \SignalNoise(t+h)}\\ \Cov{\Signal(t), X(t+h)} & = & \Cov{\Signal(t), \Signal(t+h)} \end{eqnarray}\]
We want \(S(t_0)\) from \(X(t_1), \ldots X(t_n)\) \(\Rightarrow\) we can use these covariances to find coefficients in the usual way
Wiener filter = applying these coefficients to the \(X(t)\) series to get \(\hat{\Signal}(t_0)\)
Need to know, or guess, at the correlations of either the signal or the noise to do de-noising
- Do not need any direct observations of the signal

Assume stationarity

Add the assumption of weak stationarity: \[\begin{eqnarray} \Cov{\Signal(t), \Signal(t+h)} & = & \AutoCov(h)\\ \Cov{\SignalNoise(t), \SignalNoise(t+h)} & = & \NoiseAutoCov(h)\\ \end{eqnarray}\]
Implies \(X\) is weakly stationary with \(\Cov{X(t), X(t+h)} = \AutoCov(h) + \NoiseAutoCov(h)\)
Also implies that the coefficients we use to get \(S(t_0)\) from \(X(t_1), \ldots X(t_n)\) will also work to get \(S(t_0+h)\) from \(X(t_1+h), \ldots X(t_n+h)\)
We want \(\LinearExpect{\Signal(0)| X(-r), \ldots X(0), \ldots X(r)}\)

Assume stationarity (cont’d)

Back up to just before finding the coefficients \(\vec{\beta} = (\beta_r, \ldots \beta_0, \ldots \beta_r)\): \[ \vec{\beta} \Var{\vec{X}} = \Cov{\vec{X}, S(0)} \]
Fix any \(t \in -r:r\), then \[\begin{eqnarray} \sum_{s=-r}^{r}{\beta_s \Cov{X(s), X(t)}} & = & \Cov{X(t), S(0)}\\ \sum_{s=-r}^{r}{\beta_s (\AutoCov(t-s) + \NoiseAutoCov(t-s))} & = & \AutoCov(t) \end{eqnarray}\]
As we let \(r \rightarrow \infty\), the LHS becomes convolution: \[ (\beta * (\AutoCov + \NoiseAutoCov))(t) = \AutoCov(t) \]
Fourier transforms turn convolution into multiplication: \[ \tilde{\beta}(\omega) \widetilde{(\AutoCov+\NoiseAutoCov)}(\omega) = \tilde{\AutoCov}(\omega) \]
Fourier transformation is linear: \[ \tilde{\beta}(\omega) = \frac{\tilde{\AutoCov}(\omega)}{\tilde{\AutoCov}(\omega) + \tilde{\NoiseAutoCov}(\omega)} \]

Two-sided Wiener filter for stationary processes

In summary: \[ \tilde{\beta}(\omega) = \frac{\tilde{\AutoCov}(\omega)}{\tilde{\AutoCov}(\omega) + \tilde{\NoiseAutoCov}(\omega)} = \frac{1}{1+\tilde{\NoiseAutoCov}(\omega)/\tilde{\AutoCov}(\omega)} \]
To get actual \(\beta\), inverse Fourier transform
What does this actually mean?
- Attenuate frequencies dominated by noise power
- Relatively enhance frequencies dominated by signal
\(\therefore\) not just high- or low- pass
- But measurement noise does tend to be high-frequency
  - White noise has a uniform power spectrum so it usually dominates at really high frequencies
We need to know either the noise power spectrum or the signal power spectrum
- Or their autocovariance functions (equivalent by Wiener-Khinchin theorem)
Sources:
- Experiments with known signals
- Scientific knowledge about either signal or measurement process / noise
  - E.g., 60 Hz (50 Hz in Europe)
- Guesswork
Recover \(\beta\) by inverse Fourier transform

One-sided filtering

Two-sided filtering (like we just saw) can work well for historical data (like Kyoto’s cherry blossoms)
In many practical contexts, we have \(X(-r), \ldots X(0)\) and want to guess \(S(0)\) \[ \sum_{s=-r}^{0}{\beta_s (\AutoCov(t-s) + \NoiseAutoCov(t-s))} = \AutoCov(t) \]
Or even: we have \(X(-r), \ldots X0)\) and want to guess \(S(h)\) \[ \sum_{s=-r}^{0}{\beta_s (\AutoCov(t-s) + \NoiseAutoCov(t-s))} = \AutoCov(t+h) \]
There usually isn’t a neat Fourier trick to do this exactly, but you could always invert the matrix
- Wiener worked in continuous time so he faced \(\int_{-r}^{0}{\beta(s) (\AutoCov(t-s) + \NoiseAutoCov(t-s)) ds} = \AutoCov(t+h)\), which people came to call a Wiener-Hopf integral equation

Extracting periodic or cyclic patterns

Back to the \(X(t) = \TrueRegFunc(t) + \TrueNoise(t)\) representation
Trend may be periodic, \(\TrueRegFunc(t) = \TrueRegFunc(t+\Periodicity)\) for period \(\Periodicity\)
- Also call this cyclic or seasonal
- Can also have \(\TrueRegFunc(t) = m(t) + p(t)\) with \(p\) periodic and \(m\) not, we’ll come back to that later
How can we use this, if we know \(\Periodicity\)?

Estimating the periodic component

For simplicity, say we have \(X(0), X(1), \ldots X(n-1)\), \(n=k\Periodicity\)

\[\begin{eqnarray} \EstRegFunc(0) & = & \frac{1}{k}\sum_{i=0}^{k-1}{X(i\Periodicity)}\\ \EstRegFunc(1) & = & \frac{1}{k}\sum_{i=0}^{k -1}{X(1+i\Periodicity)}\\ \EstRegFunc(t) & = & \frac{1}{k}\sum_{i=0}^{k-1}{X(t+i\Periodicity)}\\ & \vdots & \\ \EstRegFunc(\Periodicity-1) & = & \frac{1}{k}\sum_{i=0}^{k-1}{X(\Periodicity-1+i\Periodicity)} \end{eqnarray}\]

Estimating the periodic component (cont’d.)

\[\begin{eqnarray} \EstRegFunc(0) & = & \frac{1}{k}\sum_{i=0}^{k-1}{X(i\Periodicity)}\\ & = & \frac{1}{k}\sum_{i=0}^{k-1}{\TrueRegFunc(0) + \TrueNoise(i\Periodicity)}\\ & = & \TrueRegFunc(0) + \frac{1}{k}\sum_{i=0}^{k-1}{\TrueNoise(i\Periodicity)} \end{eqnarray}\]

\(\therefore\) \(\EstRegFunc(0)\) is unbiased for \(\TrueRegFunc(0)\)
Consistent if \(\Var{\EstRegFunc(0)} \rightarrow 0\), which is the same as \[ \Var{\sum_{i=0}^{k-1}{\TrueNoise(i\Periodicity)}} = \sum_{i=0}^{k-1}{\sum_{j=1}^{k-1}{\Cov{\TrueNoise(i\Periodicity), \TrueNoise(j\Periodicity)}}} = o(k^2) \]
Notice that all the covariances are for noise terms separated by at least \(\Periodicity\), if not many multiples of \(\Periodicity\), so plausibly small and shrinking

Estimating the periodic component

library(gstat)
data(wind)
dub.jan1 <- wind[wind$month == 1 & wind$day == 1, "DUB"]
mean(dub.jan1)

## [1] 12.67722

colMeans(wind[wind$month == 1 & wind$day == 1, -(1:3)])

##       RPT       VAL       ROS       KIL       SHA       BIR       DUB       CLA 
## 14.610000 11.200556 12.931667  6.845556 10.798333  7.678333 12.677222  8.952778 
##       MUL       CLO       BEL       MAL 
## 10.057222  9.936667 13.987222 18.536667

Now repeat for every day of the calendar

After extracting the periodic component

Observation = trend + fluctuation \[\begin{eqnarray} X(t) & = & \TrueRegFunc(t) + \TrueNoise(t)\\ X(t) & = & \EstRegFunc(t) + \EstNoise(t)\\ \EstNoise(t) & = & X(t) - \EstRegFunc(t) \end{eqnarray}\]

If we think \(\TrueRegFunc\) is periodic, we can estimate it by averaging over periods
Then remove the periodic component from the observations
- a.k.a. “Seasonally adjusted” values

How do we know the period \(\Periodicity\)?

Look!
Theory/prior knowledge!
Examine the power spectrum

Summary

Extracting a random signal from noise looks just like any other linear prediction problem
- \(X(t) = S(t) + \SignalNoise(t)\), predict \(S(t)\) from \(X\)
- Need to make assumptions about the noise, \(\Cov{\SignalNoise(t), \SignalNoise(t+h)}\)
Extracting periodic components by averaging
- \(X(t) = \TrueRegFunc(t) + \TrueNoise(t)\), with \(\TrueRegFunc\) deterministic
- Periodic/cyclic/seasonal component means \(\TrueRegFunc(t) = \TrueRegFunc(t+\Periodicity)\)
- Estimate the periodic trend by averaging over periods
- Need to know period \(\Periodicity\)

Time permitting: Trend + periodicity + fluctuations

\[\begin{eqnarray} X(t) & = & m(t) + p(t) + \TrueNoise(t)\\ & = & \text{long-run trend} + \text{periodic component} + \text{random fluctuations} \end{eqnarray}\]

Smooth and remove the smooth trend
- Often: deliberately under-smooth
Average over periods and remove the periodic component
What’s left is \(\approx\) fluctuations

Time permitting: Trend + Period for CO2 Concentration

Atmospheric CO2 concentration measured at the Mauna Loa observatory in Hawaii

This data set is famous as the “Keeling curve”, and goes back to the 1950s

Photo: Mauna Loa Observatory by Rick Peterson

Time permitting: Trend + Period for CO2 Concentration

Clearly some sort of repeating deviation from the long-run trend

Time permitting: Trend + Period for CO2 Concentration

Extract a linear time trend

Time permitting: Trend + Period for CO2 Concentration

Residuals from the linear trend = periodic component + fluctuations

Notice that each month clearly has its own distribution
Red line = smoothing spline applied to these data
Trough in carbon dioxide at end of northern hemisphere summer, and peak at end of northern hemisphere winter¹

Time permitting: Trend + Period for CO2 Concentration

Random fluctuations = residuals after removing the long-run trend and the periodic component

These are very structured-looking because the initial linear trend isn’t that good
- Atmospheric \(\mathrm{CO}_2\) is increasing at an accelerating rate
Iterative approach: use the initial estimate of the periodic component to get a better idea of the long-run trend

Time permitting: Trend + Period for CO2 Concentration

We could now go back to re-estimate the annual cycle, but it doesn’t change much:

Time permitting: Trend + Period for CO2 Concentration

Alternatively: if you know (or are willing to guess) a parametric form for the long-run trend, \(m(t) = f(t;\theta)\), you can introduce extra “dummy” (indicator) variables for each phase in the cycle (here, months), and use least squares to estimate \(\theta\) and the coefficients on those dummies
This is going to be somewhat more efficient, both statistically and computationally, than the iterative procedure if you have a good idea of the form of the trend
Here a reasonable guess at the long-run trend would be a quadratic, \(f(t;\theta) = \theta_0 + \theta_1 t + \theta_2 t^2\)

quad.trend.plus.cycle <- lm(co2ppm ~ date + I(as.numeric(date)^2) + factor(month),
    data = mauna_loa_monthly)
signif(coefficients(quad.trend.plus.cycle), 2)

##           (Intercept)                  date I(as.numeric(date)^2) 
##               3.3e+02               2.9e-03               1.0e-07 
##        factor(month)2        factor(month)3        factor(month)4 
##               6.4e-01               1.3e+00               2.5e+00 
##        factor(month)5        factor(month)6        factor(month)7 
##               3.0e+00               2.3e+00               5.7e-01 
##        factor(month)8        factor(month)9       factor(month)10 
##              -1.6e+00              -3.3e+00              -3.3e+00 
##       factor(month)11       factor(month)12 
##              -2.1e+00              -9.4e-01

Note: one fewer monthly dummy coefficient than the length of the cycle
- Always one fewer dummy coefficients than levels of the factor variable (do you remember why from linear regression?)
- Notice that the lowest (most negative) coefficients are for September and October, and the highest (most positive) are for April-May-June

Time permitting: Trend + Period for CO2 Concentration

The fitted values will capture the combination of the trend and the periodic cycle:

Time permitting: Trend + Period for CO2 Concentration

The residuals will be estimates of the fluctuations, which now are a lot more random looking:

Time permitting: Trend + Period for CO2 Concentration

Exercise (to try in your copious free time, not to hand in): what’d happen if you just used smooth.spline() on the data to begin with? Why might this not be a good thing to do here?

In photosynthesis, plants use carbon dioxide and water (and sunlight) to make the sugar glucose and oxygen; like us, they later metabolize that sugar for energy, which involves taking in oxygen and releasing carbon dioxide. So plants are going to remove carbon dioxide from the atmosphere, on net, when they photosynthesize more than they burn sugar. In other words, plants will be net carbon dioxide removers when they are green and growing in the spring and summer, and net carbon dioxide releasers when they aren’t, in the fall and winter. Since there are more plants in the northern hemisphere than the southern, because there’s more land mass in the northern hemisphere, the annual carbon dioxide cycle tracks the northern hemisphere seasons rather than the southern.↩

Linear Methods for Noise Removal and Filtering

In our previous episodes

Optimal linear filtering for time series

Optimal linear filtering for time series

This should look familiar

Simple case, to build intuition

One observation, one estimate

Two observations, one estimate

Two observations, one estimate

What’s going on here?

The general pattern: the Wiener filter

Assume stationarity

Assume stationarity (cont’d)

Two-sided Wiener filter for stationary processes

One-sided filtering

Extracting periodic or cyclic patterns

Estimating the periodic component

Estimating the periodic component (cont’d.)

Estimating the periodic component

After extracting the periodic component

How do we know the period \(\Periodicity\)?

Summary

Time permitting: Trend + periodicity + fluctuations

Time permitting: Trend + Period for CO2 Concentration

Time permitting: Trend + Period for CO2 Concentration

Time permitting: Trend + Period for CO2 Concentration

Time permitting: Trend + Period for CO2 Concentration

Time permitting: Trend + Period for CO2 Concentration

Time permitting: Trend + Period for CO2 Concentration

Time permitting: Trend + Period for CO2 Concentration

Time permitting: Trend + Period for CO2 Concentration

Time permitting: Trend + Period for CO2 Concentration

Time permitting: Trend + Period for CO2 Concentration

Time permitting: Trend + Period for CO2 Concentration

Time permitting: Trend + Period for CO2 Concentration