--- title: Linear Methods for Noise Removal and Filtering date: 15 September 2026 (Lecture 7) author: 36-740 output: slidy_presentation: math_method: engine: mathjax url: https://bactra.org/mathjax/tex-svg.js bibliography: locusts.bib --- ```{r, include=FALSE} # General set-up options library(knitr) opts_chunk$set(size="small", background="white", highlight=FALSE, cache=TRUE, autodep=TRUE, tidy=TRUE, tidy.opts=list(comment=FALSE), warning=FALSE, message=FALSE, echo=FALSE) ``` \[ \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 {.smaller} - 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 {.smaller} \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) {.smaller} - 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](https://www.stat.cmu.edu/~cshalizi/dst/26/lectures/04/lecture-04.html#(11)): \[ (\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 {.smaller} - 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](https://www.stat.cmu.edu/~cshalizi/dst/26/lectures/04/lecture-04.html#(15))) - 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 {.smaller} - 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.) {.smaller} \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](https://www.stat.cmu.edu/~cshalizi/dst/26/lectures/02/lecture-02.html#(29)) \[ \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 ```{r, echo=TRUE} library(gstat); data(wind) # What is the average for Dublin for January 1? dub.jan1 <- wind[wind$month==1 & wind$day==1, "DUB"] mean(dub.jan1) ## What are all the averages for January 1? colMeans(wind[wind$month==1 & wind$day==1, -(1:3)]) ``` 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 ![](https://live.staticflickr.com/4030/4456645340_56e9a569bd_z.jpg) - This data set is famous as the ["Keeling curve"](https://sioweb.ucsd.edu/programs/keelingcurve/), and [goes back to the 1950s](https://sioweb.ucsd.edu/programs/keelingcurve/2013/04/03/the-history-of-the-keeling-curve/) [Photo: Mauna Loa Observatory by Rick Peterson](https://www.flickr.com/photos/rickpeterson/4456645340/) ## Time permitting: Trend + Period for CO2 Concentration ```{r, echo=FALSE} # After http://rstudio-pubs-static.s3.amazonaws.com/237452_c385d5892dfd4e4281619346a81b4256.html mauna_loa_monthly <- read.table('ftp://aftp.cmdl.noaa.gov/products/trends/co2/co2_mm_mlo.txt') mauna_loa_monthly <- mauna_loa_monthly[, c(1, 2, 4)] names(mauna_loa_monthly) = c('year', 'month', 'co2ppm') mauna_loa_monthly$date <- as.Date(paste(mauna_loa_monthly$year, mauna_loa_monthly$month, '15', sep = '-'), format = '%Y-%m-%d') # Fix the ridiculous use of -99.99 as a missing data code mauna_loa_monthly[mauna_loa_monthly == -99.99] <- NA # Simplify life by dropping the NA rows mauna_loa_monthly <- na.omit(mauna_loa_monthly) ``` ```{r} plot(co2ppm ~ date, data=mauna_loa_monthly, pch=16, cex=0.3) ``` - Clearly some sort of repeating deviation from the long-run trend ## Time permitting: Trend + Period for CO2 Concentration - Extract a linear time trend ```{r} linear.trend <- lm(co2ppm ~ date, data=mauna_loa_monthly) plot(co2ppm ~ date, data=mauna_loa_monthly, pch=16, cex=0.3) abline(linear.trend, col="blue") ``` ## Time permitting: Trend + Period for CO2 Concentration - Residuals from the linear trend = periodic component + fluctuations ```{r} plot(residuals(linear.trend) ~ mauna_loa_monthly$month, xlab="month") annual.spline <- smooth.spline(x=mauna_loa_monthly$month, y=residuals(linear.trend)) lines(annual.spline, col="red") ``` - 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[^co2] [^co2]: 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. ## Time permitting: Trend + Period for CO2 Concentration {.smaller} - Random fluctuations = residuals after removing the long-run trend and the periodic component ```{r} plot(x=mauna_loa_monthly$date, y=residuals(annual.spline)) ``` - 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 ```{r} mauna_loa_monthly$seasonally_adjusted <- with(mauna_loa_monthly, co2ppm - predict(annual.spline, x=month)$y) plot(seasonally_adjusted ~ date, data=mauna_loa_monthly) abline(linear.trend, col="blue") smooth.trend <- with(mauna_loa_monthly, smooth.spline(x=date, y=seasonally_adjusted)) lines(smooth.trend, col="purple", lwd=2) legend("topleft", legend=c("Linear trend to raw data", "Smooth trend of adjusted data"), lty="solid", col=c("blue", "purple")) ``` ## Time permitting: Trend + Period for CO2 Concentration - We could now go back to re-estimate the annual cycle, but it doesn't change much: ```{r} mauna_loa_monthly$detrended <- with(mauna_loa_monthly, co2ppm - predict(smooth.trend, x=as.numeric(date))$y) plot(detrended ~ month, data=mauna_loa_monthly) lines(annual.spline, col="red") lines(with(mauna_loa_monthly, smooth.spline(x=month, y=detrended), col="purple")) legend("bottomleft", legend=c("Original cycle estimate (from linear trend)", "Re-estimated cycle (from smooth trend)"), lty="solid", col=c("red", "purple")) ``` ## 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$ ```{r, echo=TRUE} 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) ``` - 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: ```{r} plot(x=mauna_loa_monthly$date, y=fitted(quad.trend.plus.cycle), pch=16, cex=0.3, xlab="date", ylab="Combine quadratic trend and annual cycle") ``` ## Time permitting: Trend + Period for CO2 Concentration - The residuals will be estimates of the fluctuations, which now are a lot more random looking: ```{r} plot(x=mauna_loa_monthly$date, y=residuals(quad.trend.plus.cycle), xlab="date", ylab="Residuals from quadratic trend and annual cycle") ``` ## Time permitting: Trend + Period for CO2 Concentration ```{r} plot(annual.spline, col="red", xlab="month", ylab="annual cycle estimate", type="l") lines(with(mauna_loa_monthly, smooth.spline(x=month, y=detrended)), col="purple") lines(x=1:12, y=c(0, coefficients(quad.trend.plus.cycle)[-(1:3)]), col="green") legend("bottomleft", legend=c("Original estimate from linear trend", "Re-estimate from smooth trend", "Simulatenously estimated with quadratic trend"), lty="solid", col=c("red", "purple", "green")) ``` ## 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?