--- title: Trends and Smoothing II author: 36-467/36-667, Fall 2020 date: 10 September 2020 output: slidy_presentation --- ```{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}} \DeclareMathOperator{\tr}{tr} \DeclareMathOperator*{\argmin}{argmin} \DeclareMathOperator{\dof}{DoF} \DeclareMathOperator{\det}{det} \newcommand{\TrueNoise}{\epsilon} \newcommand{\EstNoise}{\widehat{\TrueNoise}} \] ## In our last episode... > - Data $X(t) = \TrueRegFunc(t) + \TrueNoise(t)$ > - $\TrueRegFunc$ deterministic (=**trend**), $\TrueNoise$ stochastic and mean-zero (=**fluctuations**) > - Wanted: estimates of $\TrueRegFunc$ and/or $\TrueNoise$ from _one_ data set > - Hope: $\TrueRegFunc$ is a smooth function $\Rightarrow$ average nearby $X$'s > - **Linear smoother**: $\EstRegFunc(t) = \sum_{j=1}^{n}{w(t, t_j) x_j}$ > - **Fitted values** on the data $\mathbf{\EstRegFunc} = \mathbf{w}\mathbf{x}$ > - $\mathbf{w}$ determines the properties of the fitted values > + $\mathbf{w}$ has $n$ eigenvalues $\lambda_1 \geq \lambda_2 \geq \ldots \geq \lambda_n$, and eigenvectors $\mathbf{v}_1, \mathbf{v}_2, \ldots \mathbf{v}_n$ > + Components of $\mathbf{x}$ which look like eigenvectors $\mathbf{v}$ with big eigenvalues $\lambda$ will be preserved by the smoothing ($\lambda = 1$), or only shrunk a little ($\lambda$ nearly 1) > + Components of $\mathbf{x}$ which look like eigenvectors with small eigenvalues will be shrunk a lot by smoothing ($\lambda$ near zero) or completely eliminated ($\lambda = 0$) > + The smoother is biased towards finding trends that look like its leading eigenvectors > + degrees of freedom $=$ covariance between data and fitted values $=$ sum of the eigenvalues of $\mathbf{w}$ ## An example to make things concrete: GDP Per Capita - GDP = gross domestic product = annual total value money value of all goods and services sold in a country - 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)) 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: 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="blue") ``` - 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 ## 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 > - Every entry in $\mathbf{w}$ is $1/k$ (here $1/5$) or 0 > + $1/k$ along the diagonal (why?) > - 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) > - After that we get sine waves (as we saw last time) > + Bigger eigenvalues $\Leftrightarrow$ longer wavelength sine waves > + "low-pass filter" > - The moving average "likes" (preserves, is biased towards) patterns in the data that look like 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="purple", lwd=2) legend("bottomleft", legend=c("Data", "MA(5)", "MA(21)"), lty="solid", col=c("black", "red", "purple")) ``` ## 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="darkgreen", 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") ``` ## 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") ``` ## 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") ``` ## 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} > - _Convince yourself_: $\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 sometimes called 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 more like the price at which we'll trade more curvature ($m^{\prime\prime}$) for less mean squared error > - 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) ## 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} ## Many variants - $h$-block CV: omit a buffer of radius $h$ around the hold-out point from the training set - $k$- or $v$-fold CV: divide data into $k$ equal-sized "folds", try to predict each fold using the rest of the data - $hv$-block CV: $v$-fold with a buffer - etc., etc. ## The moral on CV - _Never_ care about how good the in-sample fit is ($R^2$, $R^2_{adj}$, etc.) - _Always_ care about ability to predict new data ## Spline smoothing of economic growth - `smooth.spline()` is pretty robust but chokes on NA values: ```{r, echo=TRUE} growth.ss <- with(na.omit(gdppc), 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="purple", lwd=2) ``` ## 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 to smooth by cross-validation