--- title: Trends and Smoothing date: 27 August 2026 (Lecture 2) author: "36-740, Fall 2026" 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, warning=FALSE, message=FALSE, echo=FALSE) ``` \[ \newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]} \newcommand{\Var}[1]{\mathrm{Var}\left[ #1 \right]} \newcommand{\Cov}[1]{\mathrm{Cov}\left[ #1 \right]} \newcommand{\TrueRegFunc}{\mu} \newcommand{\EstRegFunc}{\widehat{\TrueRegFunc}} \newcommand{\TrueNoise}{\epsilon} \newcommand{\EstNoise}{\widehat{\TrueNoise}} \DeclareMathOperator{\tr}{tr} \DeclareMathOperator*{\argmin}{argmin} \DeclareMathOperator{\dof}{DoF} \] ## The mathematical setting - Treat the data as the realization of a **random process** or **stochastic process** = random variables spread over time and space - At point $r$ and time $t$, we see $X(r,t)$ + Just $X(r)$ if only spatial, just $X(t)$ if only temporal - Generally, $X(r,t)$ and $X(q,s)$ are statistically dependent - Can also think of the process as a _random function_ of the coordinates ## Trend = central tendency of the random process \[ \TrueRegFunc(r,t) \equiv \Expect{X(r,t)} \] - This is what we'll mean by the **trend** - $\TrueRegFunc$ is a _deterministic_ function of the coordinates; not random ## Trend + Fluctuations Random process = trend + fluctuations \begin{eqnarray} X(r,t) & = & \TrueRegFunc(r,t) + \TrueNoise(r,t)\\ \Expect{\TrueNoise(r,t)} & = & 0 \end{eqnarray} (Sometimes called "signal plus noise" representation) $\TrueNoise =$ another random process - Can't, _in general_, assume more about $\TrueNoise$: + It might be correlated + It might not have constant variance + It might have different distributions everywhere + It might not be Gaussian + etc., etc. - But we _can_ assume $\Expect{\TrueNoise(r,t)} = 0$ ## Learning the Trend - We want to learn $\TrueRegFunc$ - Either we have many realizations of the process or we just have one - Things are a lot easier if we have many realizations ## Averaging independent realizations gives the trend - We run the experiment many times, getting independent realizations $X^{(1)}, X^{(2)}, \ldots X^{(m)}$ - Extracting the trend is simple, because the law of large numbers applies: \[ \frac{1}{m}\sum_{i=1}^{m}{X^{(i)}(r,t)} \rightarrow \TrueRegFunc(r,t) \] - Used a lot in neuroscience, physics, chemistry... - Harder to use when we're not doing controlled experiments ## Trends from a _single_ realization - We only have one history of cherry-blossoms in Kyoto + Or US economy post-WWII or shale deposits under Pennsylvania, etc. - To get a trend, we'll need to assume _something_ - Assume a whole theory - _or_ Assume a little about $\TrueRegFunc$ ## Assume a whole theory {.smaller} - Sometimes we actually have good theories that tell us what the trend should be \[ \TrueRegFunc(r, t) = f(r,t;\theta) \] with $f$ completely specified except for some constants $\theta$ - Then all we have to do is estimate $\theta$ - Usual approach is least squares: \[ \hat{\theta} = \argmin_{\theta}{\left(\frac{1}{n}\sum_{i=1}^{n}{(X(r_i, t_i) - f(r_i, t_i;\theta))^2} \right)} \] + Issues with correlations in $\TrueNoise$ that we'll come back to later in the course - Situations where this applies: + Physics, astronomy * This is where the method of least squares came from [@Farebrother-on-linear-relationships] + Some other natural sciences: chemistry, geology, climatology... + Parts of biology: population genetics, ecology, epidemiology + Social sciences: demography, parts of economics (maybe) ## Assume a whole theory - Pros: + We know how to estimate parameters! + We can interpret estimated parameters using the theory + Statistically efficient if the model is right + Good for extrapolating if the model is right - Cons: + Doesn't check the theory + Systematically misleading if the model is wrong + Only works if the theory is "tight" enough to give us a function up to unknown constants; most of science isn't like that ## Assume something about $\TrueRegFunc$ > - Assume $\TrueRegFunc$ is a **smooth** function > - Concretely: assume $\TrueRegFunc$ changes fairly slowly, $\TrueRegFunc(t-h) \approx \TrueRegFunc(t) \approx \TrueRegFunc(t+h)$ > - We know that $X(t)$ is a noisy observation of $\TrueRegFunc(t)$ > - With this assumption, $X(t-h)$ and $X(t+h)$ are _also_ noisy observations of $\TrueRegFunc(t)$ (plus a little bit more error) > - $\therefore$ If $\TrueRegFunc$ is smooth, we can estimate it by averaging nearby observations ## Our first smoothing estimator \[ \EstRegFunc(t_i) = \frac{1}{3}\sum_{k=-1}^{k=+1}{X(t_{i+k})} \] - Why do this? \begin{eqnarray} \EstRegFunc(t_i) & = & \frac{1}{3}\sum_{k=-1}^{k=+1}{\TrueRegFunc(t_{i+k})} + \frac{1}{3}\sum_{k=-1}^{k=+1}{\TrueNoise(t_{i+k})}\\ & = & \TrueRegFunc(t_i) + \frac{1}{3}\sum_{k=-1}^{k=+1}{(\TrueRegFunc(t_{i+k}) - \TrueRegFunc(t_i))} + \frac{1}{3}\sum_{k=-1}^{k=+1}{\TrueNoise(t_{i+k})}\\ & = & \text{truth} + \text{bias} + \text{noise} \end{eqnarray} - Bias term: + Big if $\TrueRegFunc$ changes rapidly and observations are far apart + Small if $\TrueRegFunc$ changes slowly and observations are close - Noise: mean 0, usually smaller than the noise variance in any one observation ## Linear smoothers - Data: $x(t_1), \ldots x(t_n)$; abbreviate as $x_1, \ldots x_n$. - A **linear smoother** is an estimator that looks like this: \[ \EstRegFunc(t) = \sum_{j=1}^{n}{w(t, t_j) x_j} \] - Notice: + Linear in the observations + Can be _very_ nonlinear in $t$ + Same idea with space but more notation - Different linear smoothers $=$ different choices of weights $w$ ## Examples of linear smoothers {.smaller} - _Global constant_: The trivial case, $w=1/n$ no matter what. - $k$ _Nearest neighbors_ (kNN): $w(t, t_j) = 1/k$ if $t_j$ is one of the $k$ points closest to $t$, and otherwise $w=0$. - _Moving average_ (MA): $w(t, t_j) = 0$ outside some distance from $t$, and constant within it. + Also one-sided moving averages (usually only from the past) - _Weighted moving averages_: $w(t, t_j) = g(|t-t_j|/\tau)$, for some fixed function $g$ + Often $g$ is the (negative) exponential function, so **exponentially-weighted moving averages**. - _Kernels_: like a weighted moving average, but weights are **normalized** (non-negative, always sum to 1), like picking a pdf $K$ and a **bandwidth** $h$ and setting \[ w(t, t_j) = \frac{K(|t-t_j|/h)}{\sum_{j^{\prime}=1}^{n}{K(|t-t_{j^{\prime}}|/h)}} \] + "Attention" in transformers/large language models [@Attention-is-all-you-need] is one-sided kernel smoothing with a particular, even peculiar, kernel [@Tsai-et-al-attention-is-kernel-smoothing] - _Linear and polynomial regression on the coordinates_: Linear smoothers but with very weird weights - _Splines_: Smooth piecewise polynomials that actually do what you hoped polynomial regression would - _Neural networks_: Output of the last layer is linear in features defined by the lower layers ## Properties of linear smoothers - **Fitted values**: at each $t_i$ where we got data, we have $\EstRegFunc(t_i)$, for short $\EstRegFunc_i$ - Group the fitted values in to a vector ($n\times 1$ matrix): \[ \mathbf{\EstRegFunc} = \left[\begin{array}{c} \EstRegFunc(t_1) \\ \EstRegFunc(t_2) \\ \vdots \\ \EstRegFunc(t_n) \end{array} \right] \] - Group all the observations of $x$ into a vector $\mathbf{x}$ - The fitted values are the observations _times_ a matrix: \[ \mathbf{\EstRegFunc} = \mathbf{w} \mathbf{x} \] where $w_{ij} = w(t_i, t_j)$. - The matrix $\mathbf{w}$ is the **smoothing**, **influence** or **hat** matrix - Important: the fitted values $\mathbf{\EstRegFunc}$ are _random_ ## Expectation of the fitted values; bias - Expected fitted values: \begin{eqnarray} \Expect{\mathbf{\EstRegFunc}} & = & \Expect{\mathbf{w}\mathbf{X}}\\ & = & \mathbf{w}\Expect{\mathbf{\TrueRegFunc} + \mathbf{\TrueNoise}}\\ & = & \mathbf{w}\Expect{\mathbf{\TrueRegFunc}} + \mathbf{w}\Expect{\mathbf{\TrueNoise}}\\ & = & \mathbf{w}\mathbf{\TrueRegFunc} \end{eqnarray} - When is a linear smoother unbiased? \begin{equation} \text{unbiased trend estimates} \Leftrightarrow \mathbf{\TrueRegFunc} = \mathbf{w}\mathbf{\TrueRegFunc} \end{equation} - A linear smoother is unbiased if, and only if, $\mathbf{\TrueRegFunc}$ is an eigenvector of $\mathbf{w}$ with eigenvalue 1 ## Shrinkage {.smaller} - For (most) matrices, the eigenvectors form a basis: for any vector $\mathbf{x}$, there are $n$ coefficients $c_i$ s.t. \[ \mathbf{x} = \sum_{i=1}^{n}{c_i \mathbf{v}_i} \] + We break the data up in to different "components" defined by the influence matrix $\mathbf{w}$ - After smoothing: \[ \mathbf{w}\mathbf{x} = \sum_{i=1}^{n}{\mathbf{w} c_i \mathbf{v}_i} = \sum_{i=1}^{n}{c_i \lambda_i \mathbf{v}_i} \] - Components of the data which match eigenvectors $\mathbf{v}_i$ with big eigenvalues $\lambda_i$ get enhanced/amplified by smoothing - Components of the data which match eigenvectors with small eigenvalues get **shrunk** + If $\lambda_i = 0$, that component of the data is completely eliminated by smoothing - Usually, the biggest eigenvalue of $\mathbf{w}$ will be 1 + For a linear trend (=linear regression on $t$), every eigenvalue is either 0 or 1 + For (weighted) averaging, the largest eigenvalue will be 1 ## A little example {.smaller} - Smoothing by averaging to either side + Need to handle the end-points... ```{r, echo=TRUE} n <- 10 w <- matrix(0, nrow=10, ncol=10) diag(w) <- 1/3 for (i in 2:(n-1)) { w[i,i+1] <- 1/3 w[i,i-1] <- 1/3 } w[1,1] <- 1/2 w[1,2] <- 1/2 w[n,n-1] <- 1/2 w[n,n] <- 1/2 ``` ## A little example {.smaller} - Eigenvalues: ```{r} eigen(w)$values ``` - Leading eigenvector: ```{r} eigen(w)$vectors[,1] ``` - This is _constant_, $\propto \mathbf{1}$ + If every entry in $\mathbf{x}$ is the same, averaging changes nothing! ## Eigenvectors 2--4 of the example ```{r, echo=FALSE} plot(1:10, eigen(w)$vectors[,2], lty="solid", type="l", xlab="Position", ylab="Eigenvectors of smoothing matrix") lines(1:10, eigen(w)$vectors[,3], lty="dashed") lines(1:10, eigen(w)$vectors[,4], lty="dotted") rounded.eigenvalues <- signif(eigen(w)$values[2:4], 2) legend("bottomleft", legend=c(paste("lambda=", rounded.eigenvalues[1]), paste("lambda=", rounded.eigenvalues[2]), paste("lambda=", rounded.eigenvalues[3])), lty=c("solid","dashed","dotted")) ``` ## Variance of the fitted values - Remember, fitted values are random: \begin{eqnarray} \Var{\mathbf{\EstRegFunc}} & = & \Var{\mathbf{w}\mathbf{X}}\\ & = & \mathbf{w}\Var{\mathbf{X}}\mathbf{w}^T\\ & = & \mathbf{w}\Var{\mathbf{\TrueRegFunc} + \mathbf{\TrueNoise}}\mathbf{w}^T\\ & = & \mathbf{w}\Var{\mathbf{\TrueNoise}}\mathbf{w}^T \end{eqnarray} - Unfortunately we don't know $\Var{\mathbf{\TrueNoise}}$ - **IF** we can assume $\Var{\mathbf{\TrueNoise}} = \sigma^2 \mathbf{I}$, then \[ \Var{\EstRegFunc} = \sigma^2 \mathbf{w}\mathbf{w}^T \] as a special case ## Degrees of freedom {.smaller} - How much do the fitted values co-vary with the data? \begin{eqnarray} \sum_{i=1}^{n}{\Cov{\EstRegFunc_i, X_i}} & = & \sum_{i=1}^{n}{\Cov{\sum_{j=1}^{n}{w_{ij} X_j}, X_i}}\\ & = & \sum_{i=1}^{n}{\sum_{j=1}^{n}{w_{ij} \Cov{X_i, X_j}}}\\ & = & \sum_{i=1}^{n}{\sum_{j=1}^{n}{w_{ij} \Cov{\TrueNoise_i, \TrueNoise_j}}} \end{eqnarray} - **IF** $\Cov{\TrueNoise_i, \TrueNoise_j} = 0$ if $i\neq j$, and $=\sigma^2$ if $i=j$, _then_ as a special case \[ \sum_{i=1}^{n}{\Cov{\EstRegFunc_i, X_i}} = \sigma^2 \tr{\mathbf{w}} \] - People use this to _define_ the **effective number of degrees of freedom**: \[ \dof{\mathbf{w}} \equiv \tr{\mathbf{w}} \] + Remember $\tr{\mathbf{w}} \equiv \sum_{i=1}^{n}{w_{ii}}$ + Also remember $\tr{\mathbf{w}} =$ sum of the eigenvalues of $\mathbf{w}$ ## Using linear smoothers - Pros: + Don't need a theory about what the trend should be + Can _discover_ the underlying trend + Can be used to _check_ theories which make predictions about the trend + Mostly straightforward calculations - Cons: + Doesn't use the data as efficiently as a theory-based model * ...when the model is right + Need to decide _how much_ smoothing/averaging to do + Need to decide _what kind_ of smoothing/averaging to do * Usually less important than the "how much" question + Extrapolation is hazardous ## An example to make things concrete: GDP Per Capita - GDP = gross domestic product = total money value of all goods and services sold in a country per year - GDP per capita = average annual income per person - US GDP per capita, adjusted for inflation, measured quarterly: ```{r} library(pdfetch) # Fetch real (inflation-adjusted) US per-capita gross domestic product # ("chained 2012 dollars", i.e., no inflation adjustment for 2012) gdppc.fred <- pdfetch_FRED("A939RX0Q048SBEA") # This comes as a complicated data type, so break it down to a simple data # frame, converting dates to a year with a decimal fraction library(xts) # Functions needed for pdfetch's preferred format library(lubridate) # Provides useful date-conversion functions gdppc <- data.frame(year=decimal_date(index(gdppc.fred)), y=as.numeric(gdppc.fred)) ``` ```{r, fig.gdp.logscale} plot(gdppc, type="l", ylab="GDPPC (dollars/person/yr)", main="Real GDP per capita (2012 dollars)", ylim=c(0, max(gdppc))) ``` ## An example to make things concrete: GDP Per Capita - GDP = gross domestic product = total money value of all goods and services sold in a country per year - GDP per capita = average annual income per person - US GDP per capita, adjusted for inflation, measured quarterly: ```{r, fig.gdp} plot(gdppc, type="l", ylab="GDPPC (dollars/person/yr)", main="Real GDP per capita (2012 dollars)", ylim=c(1e4, max(gdppc)), log="y") ``` ## An example to make things concrete: GDPPC growth rate - Logarithmic growth rate: ```{r} gdppc$growth <- c(NA, diff(log(gdppc$y))*4) ``` ```{r, fig.gdp.growth} plot(growth ~ year, data=na.omit(gdppc), type="l", xlab="year", ylab="Growth rate (1/yr)") ``` ## Trend=mean? - One possible trend = constant at the global mean + Interpretation: steady growth plus random fluctuations ```{r} <> mean.growth <- mean(gdppc$growth, na.rm=TRUE) abline(h=mean.growth, col="green") ``` - Might make more sense to omit some extreme values when estimating the mean? ## The math for taking the trend to be the global mean > - $\mathbf{w}$ is the $n\times n$ matrix with $1/n$ everywhere > - You can check: eigenvalues are 1 and 0 (repeated $n-1$ times) > - You can check: eigenvector for $\lambda=1$ is $[1 1 \ldots 1]$ > - Degrees of freedom = 1 ## A Law of Large Numbers > - Suppose $\mu_t=\mu$ for all $t$ (no trend, constant trend) > - Look at $A_n = n^{-1}\sum_{t=1}^{n}{X_t}$, the **time average** > - $\Expect{(A_n - \mu)^2} = (\Expect{A_n} - \mu)^2 + \Var{A_n}$ > - $\Expect{A_n} = \frac{n\mu}{n} = \mu$, an unbiased estimate > - $\Var{A_n} = n^{-2}\left(\sum_{t=1}^{n}{\Var{X_t}} + \sum_{t=1}^{n}\sum_{s\neq t}{\Cov{X_t, X_s}}\right) \equiv V_n/n^2$ > - $V_n = o(n^2)$ $\Leftrightarrow$ $\Expect{(A_n - \mu)^2} \rightarrow 0$ > - If $\Var{X_t}$ does not blow up, and if $\Cov{X_t, X_{t+h}} \rightarrow 0$ as $|h| \rightarrow \infty$ fast enough, for all $t$, then $V_n = o(n^2)$ > - Goes back to @Taylor-diffusion-by-continuous-movements, see @CRS-simplest-ergodic-theorem ## Trend=short moving average? - Moving average going out 0.5 yr (=2 quarters) on either side: ```{r} <> library(zoo) ma.5 <- rollmean(gdppc$growth, k=(2*2+1), fill=NA) lines(gdppc$year, ma.5, col="red", lwd=2) ``` ## The math for moving averages {.smaller} > - Every entry in $\mathbf{w}$ is $1/k$ (here $1/5$) or 0 > - Degrees of freedom = $\tr{\mathbf{w}} = n/k$ > - Leading eigenvalue is always 1 > + Generally true when we do (weighted) averages; see handout > + Other eigenvalues get smaller as width of the averaging window grows > - Leading eigenvector is constant (like with the global mean), $\mathbf{1}$ > - After that we get sine waves > + Bigger eigenvalues $\Leftrightarrow$ longer wavelength sine waves > + "low-pass filter" > - The moving average "likes" (preserves, is biased towards) patterns in the data that resemble slowly-changing sine waves > + Or sums of slowly-changing sine waves > + This bias gets stronger as the width of the averaging grows ## Trend=wider moving average? - Moving average going out 2.5 yr (=10 quarters) on either side + Correlations are pretty small after $\approx 10$ quarters ```{r} <> ma.21 <- rollmean(gdppc$growth, k=(10*2+1), fill=NA) lines(gdppc$year, ma.5, col="red", lwd=2) lines(gdppc$year, ma.21, col="orange", lwd=2) legend("bottomleft", legend=c("Data", "MA(5)", "MA(21)"), lty="solid", col=c("black", "red", "orange")) ``` ## Trend=one-sided moving average? - Might seem weird to have 2020 affecting trend for 2018... - Moving average over the previous four years (=16 quarters) ```{r} <> ma.16.past <- rollmean(gdppc$growth, k=16, fill=NA, align="right") lines(gdppc$year, ma.16.past, col="magenta", lwd=2) ``` ## Data = trend + fluctuation > - $X(t) = \TrueRegFunc(t) + \TrueNoise(t)$ > - $\Rightarrow$ $\TrueNoise(t) = X(t) - \TrueRegFunc(t)$ > - $\Rightarrow$ $\EstNoise(t) \equiv X(t) - \EstRegFunc(t) =$ **residuals** ## Some residuals Residuals from using a constant trend: ```{r} plot(gdppc$year, gdppc$growth-mean.growth, xlab="year", ylab="fluctuation in growth rate (1/yr)", type="l", col="green") ``` ## Some residuals Residuals from using an MA(5): ```{r} plot(gdppc$year, gdppc$growth-ma.5, xlab="year", ylab="fluctuation in growth rate (1/yr)", type="l", col="red") ``` ## Some residuals Residuals from using an MA(16) on the past: ```{r} plot(gdppc$year, gdppc$growth-ma.16.past, xlab="year", ylab="fluctuation in growth rate (1/yr)", type="l", col="magenta") ``` ## Some math of residuals \begin{eqnarray} \mathbf{\EstNoise} & = & \mathbf{x} - \mathbf{\EstRegFunc}\\ & = & \mathbf{x} - \mathbf{w}\mathbf{x}\\ & = & (\mathbf{I} - \mathbf{w})\mathbf{x} \end{eqnarray} $\mathbf{I}-\mathbf{w}$ has same eigenvectors as $\mathbf{w}$, but eigenvalues $1-\lambda$ ## Expected residuals \begin{eqnarray} \Expect{\mathbf{\EstNoise}} & = & \Expect{(\mathbf{I}-\mathbf{w})\mathbf{X}}\\ & = & (\mathbf{I}-\mathbf{w})\mathbf{\TrueRegFunc} \end{eqnarray} Biased trend estimate $\Leftrightarrow$ biased fluctuation estimate ## Variance and covariance of the residuals \[ \Var{\mathbf{\EstNoise}} = (\mathbf{I}-\mathbf{w}) \Var{\mathbf{\epsilon}} (\mathbf{I}-\mathbf{w})^T \] IF $\Var{\mathbf{\epsilon}} = \sigma^2 \mathbf{I}$, THEN $\Var{\mathbf{\EstNoise}}= \sigma^2 (\mathbf{I}-\mathbf{w})(\mathbf{I}-\mathbf{w})^T$ **NB**: Correlations from off-diagonal entries in $\mathbf{w}$, _even though_ there are no correlations for the true fluctuations - The way that smoothing creates correlations in fitted and detrended values is the **Yule-Slutsky effect** + See the handout for the origin of the name and more details ## Splines \[ \EstRegFunc = \argmin_{m}{\frac{1}{n}\sum_{i=1}^{n}{(x_i - m(t_i))^2} + \lambda\int{(m^{\prime\prime}(t))^2 dt}} \] > - This $\lambda$ not an eigenvalue (sorry) > + it's the price at which we'll trade more curvature ($m^{\prime\prime}$) for less mean squared error > + equivalently (Lagrange multipliers), minimize MSE subject to a constraint on the curvature > + equivalently (duality), minimize curvature subject to a constraint on the MSE > - Fit to the data points vs. over-all curvature > - Minimization is over _all_ functions > - Solution is always a piecewise cubic polynomial, but continuous, with continuous 1st and 2nd derivatives > - $\lambda \rightarrow 0$ $\Rightarrow$ Straight lines between data points > - $\lambda \rightarrow \infty$ $\Rightarrow$ Global linear fit > - $\downarrow$ degrees of freedom as $\uparrow \lambda$ > - Easiest R command: `smooth.spline()` (described in detail in the handout) > + Many, many more advanced packages... ## How do we pick $\lambda$? > - Want trend to predict not-yet-seen stuff (interpolate, extrapolate, filter) > - A good $\lambda$ predicts new stuff well > - Hold out part of the data and try to predict that from the rest ## Leave-one-out cross-validation (LOOCV) - For each of the $n$ data points: + Fit using every data point _except_ $i$, get $\EstRegFunc^{(-i)}$; + Find $\EstRegFunc^{(-i)}(t_i)$; + Find $(x_i - \EstRegFunc^{(-i)}(t_i))^2$. - Average over all data points, $n^{-1}\sum_{i=1}^{n}{(x_i - \EstRegFunc^{(-i)}(t_i))^2}$ - Low LOOCV $\Leftrightarrow$ good ability to predict new data - This is what `smooth.spline` does automatically ## Leave-one-out cross-validation (LOOCV) Don't have to re-fit linear smoothers $n$ times: \begin{eqnarray} \EstRegFunc^{(-i)}(t_i) &= & \frac{({\mathbf{w} \mathbf{x})}_i - w_{ii} x_i}{1-w_{ii}}\\ x_i - \EstRegFunc^{(-i)}(t_i) & = & \frac{x_i - \EstRegFunc(t_i)}{1-w_{ii}}\\ LOOCV & = & \frac{1}{n}\sum_{i=1}^{n}{\left(\frac{x_i-\EstRegFunc(t_i)}{1-w_{ii}}\right)^2} \end{eqnarray} This just needs _one_ fit to the whole data, and the diagonal of $\mathbf{w}$ (Trick due to Grace Wahba in the 1970s, see @Wahba-spline-models, Thm. 4.2.1) ## Spline smoothing of economic growth ```{r, echo=FALSE} # smooth.spline() is pretty robust but chokes on NA values # Let's omit the weirdness of 2020--2021 growth.ss <- with(na.omit(gdppc[(gdppc$year < 2020) | (gdppc$year > 2021),]), smooth.spline(x=year, y=growth)) growth.ss ``` ```{r} <> # The output of smooth.spline() includes $x and $y components designed to # work well with plotting commands line points() and lines() lines(growth.ss, col="blue", lwd=2) ``` ## Spline smoothing of economic growth Confession: I cheated a little by omitting 2020--2021 ```{r, echo=FALSE} # Let's include the weirdness growth.ss.all <- with(na.omit(gdppc), smooth.spline(x=year, y=growth)) growth.ss.all ``` ```{r} <> lines(growth.ss, col="blue", lwd=2) lines(growth.ss.all, col="cyan", lwd=2) legend("bottomleft", legend=c("Omitting 2020--1", "Including 2020--1"), lty="solid", pch=NA, col=c("blue", "cyan")) ``` - Q: Could we omit 2020--2021 from the cross-validation, but still include it in the estimation? - A: Yes. (_Not shown in class_) Here's the result: ```{r} growth.ss.hybrid <- with(na.omit(gdppc), smooth.spline(x=year, y=growth, lambda=growth.ss$lambda)) ``` ```{r} <> lines(growth.ss, col="blue", lwd=2) lines(growth.ss.all, col="cyan", lwd=2) lines(growth.ss.hybrid, col="cyan4", lwd=2) ``` Unsurprisingly (?), this looks almost the same as the original curve, _except_ around 2020 ## Summing up > - If the trend is smooth, we can estimate it by smoothing > - Every smoother is biased towards some patterns and against others > - Properties of the fitted values come from the weights > - Fluctuations are estimated as residuals after removing a trend > - De-trending can create correlations > - We decide how much to smooth by cross-validation (usually) ## Next time > - Breaking down into components that are _not_ defined by our smoothing method > - Getting smoothing anyway by strategic choice of which components to shrink ## Exercises, to think through / practice with, rather than hand in {.smaller} 1. Extract residuals from splines (i.e., both the spline curves) and plot them. Qualitatively, how do they look compared to each other? Compared to the residuals from the moving averages? 2. Make a mosaic of plots of the residuals from the global mean, from all the moving averages, and from the splines, against each other. Qualitatively, do they mostly line up around the diagonal, or are there visually-discernable patterns? + Make sure you're plotting residuals for the same time points across different estimators; this may need some fiddly adjustment of indices 3. ([_Spiders Georg_](https://www.tumblr.com/reallyreallyreallytrying/40033025233/average-person-eats-3-spiders-a-year-factoid)) Give 1--3 reasons why it makes sense to not include 2020--2021 in cross-validation here ("adn should not have been counted"). Also give 1--3 reasons why it is important, even vital, to include those years. 4. Refer to the slide "A Law of Large Numbers". Suppose $\Expect{X_t} = \mu$, $\Var{X_t} = \sigma^2$, and $\Cov{X_t, X_s} = \sigma^2 \beta^{|t-s|}$ for some $\beta \in (-1, 1)$. Show that $V_n = o(n^2)$ and find an expression for $\lim{n \Expect{(A_n - \mu)^2}}$ in terms of $\mu$, $\sigma^2$ and/or $\beta$. + See also the exercises at the end of lecture 1. 5. (LLN continued) More generally, suppose that $\Cov{X_t, X_s} = \rho(|t-s|)$, with $\rho(0) = \sigma^2$. Show that if $\sum_{h=-\infty}^{\infty}{\rho(h)} = \sigma^2 \tau < \infty$, then $V_n/n \rightarrow \sigma^2 \tau$ and $n\Expect{(A_n - \mu)^2} \rightarrow \sigma^2 \tau$. Explain why, in this situation $n/\tau$ is called the **effective sample size**. ## Backup: More about eigenvalues and eigenvectors of influence matrices {.smaller} ### Linear Regression $\Leftrightarrow$ Projection onto a linear subspace - Say we only have one coordinate $t$, so data point $i$ was observed at $t_i$. Collect all the $t_i$ into a vector ($n\times 1$ matrix)$ $\mathbf{t}$. Then in linear regression, the fitted values _must_ take the form \[ \mathbf{\EstRegFunc} = \beta_0 \mathbf{1} + \beta_1 \mathbf{t} \] - $\therefore$ Least squares takes the data vector $\mathbf{x}$ and projects into the linear subspace spanned by $\mathbf{1}$ and $\mathbf{t}$ - If (miraculously) $\mathbf{x}$ is already in that subspace, _nothing changes_ - $\therefore$ Projection onto a linear subspace is an **idempotent** linear operator, where $\mathbf{w}\mathbf{w} = \mathbf{w}$ $\Rightarrow$ the only possible eigenvalues are 0 and 1 - Say there are $d$ coordinates, $r^{(1)}, \ldots r^{(d)}$. Form the vectors $\mathbf{r}^{(1)}, \ldots \mathbf{r}^{(d)}$. Linear regression still enforces \[ \mathbf{\EstRegFunc} = \beta_0 \mathbf{1} + \beta_1 \mathbf{r}^{(1)} + \ldots \beta_d \mathbf{r}^{(d)} \] so we're still projecting $\mathbf{x}$ into a $(d+1)$-dimensional linear subspace, and the only possible eigenvalues are still 1 and 0 - Say we linearly regress on $k$ (nonlinear) transformations of the coordinates, $f_1(\mathbf{r}, t), \ldots f_k(\mathbf{r}, t)$. Form the vectors $\mathbf{f}^{(1)}, \ldots \mathbf{f}^{(k)}$ and apply the same logic. ## Backup: More about eigenvalues and eigenvectors of influence matrices {.smaller} #### Weighted averages - Suppose the linear smoother always predicts a weighted average of the data values, so $\sum_{j=1}^{n}{w(t, t_j)} = 1$, $w(t, t_j) \geq 0$ for all $t$. - It follows that $0 \leq w_{ij} \leq 1$ for any such smoother. (For each $i$, if $w_{ij}=1$ for some $j$, $w_{ik}=0$ for all $k\neq j$.) - If $x_i = x$ for all $i$, then $\EstRegFunc_i = \sum_{j=1}^{n}{w(t_i, w_j) x} = x\sum_{j=1}^{n}{w(t_i, t_j)} = x$. Therefore $\mathbf{1}$ must be an eigenvector of $\mathbf{w}$ with eigenvalue $1$. + Because eigenvectors with distinct eigenvalues are orthogonal, all the other eigenvectors must have both positive and negative entries. - To prove there's no eigenvalue $\lambda >1$, suppose there was one, with eigenvector $\mathbf{v}$. We'd have $\lambda v_i = \sum_{j}{w_{ij} v_j}$ for all $i$. Now say that $v_{\min} < 0$ and $v_{\max} > 0$ are the smallest and largest entries in $\mathbf{v}$, at positions $m$ and $M$. Because $0 \leq w_{ij} \leq 1$, $\sum_{j}{w_{ij} v_j}$ must lie between $v_{\min}$ and $v_{\max}$ for all $i$. But now we have a contradiction: $\lambda v_{\max} > v_{\max}$, because $\lambda > 1$, but also $\lambda v_{\max} = \sum_{j}{w_{M j} v_j} \leq v_{\max}$. + Similarly, but with $v_{\min}$, for showing that there can't be an eigenvalue $< -1$. - Repeated _averaging_ shrinks towards the subspace of the constant vector $\mathbf{1}$, with different directions shrunk at different rates + "Softer" than immediately projecting into a linear subspace (which contains $\mathbf{1}$) ## References