How can we explain the "bad reputation" of higher-order polynomials? We all must have heard it by now - when we start learning about statistical models overfitting data, the first example we are often given is about "polynomial functions" (e.g., see the picture here):

We are warned that although higher-degree polynomials can fit training data quite well, they surely will overfit and generalize poorly to the test data.
Why does this happen? Is there a mathematical justification as to why (higher-degree) polynomial functions overfit the data? The closest explanation I could find online was something called "Runge's phenomenon", which suggests that higher-order polynomials tend to "oscillate" a lot - does this explain why polynomial functions are known to overfit data?
I understand that there is a whole field of "regularization" that tries to fix these overfitting problems (e.g., penalization can prevent a statistical model from "hugging" the data too closely) - but just using mathematical intuition, why are polynomials known to overfit the data?
In general, "functions" (e.g., the response variable you are trying to predict using machine learning algorithms) can be approximated using older methods like Fourier series, Taylor series and newer methods like neural networks. I believe that there are theorems that guarantee that Taylor series, polynomials and neural networks can "arbitrarily approximate" any function. Perhaps neural networks can promise smaller errors for simpler complexity?
But are there mathematical reasons behind higher-order polynomials (e.g., polynomial regression) being said to have a bad habit of overfitting, to the extent that they have become very unpopular? Is this solely explainable by Runge's phenomenon?
Reference:
Gelman, A. and Imbens, G. (2019) Why high order polynomials should not be used in regression discontinuity designs. Journal of Business and Economic Statistics 37(3), pp. 447-456. (An NBER working paper version is available here)

You have more problems than just Runge's phenomenon. Below is an example for fitting a tenth degree polynomial to  21 data points that follow the curve
$$y = sin(6\pi x^2)$$

*

*Runge's phenomenon: The black broken line is the least-squares fit to these 21 points if there is no noise. You see some larger error towards the edges, but when you only interpolate then it is not incredibly large.


*Variance and overfitting. On top of that phenomenon, you get that when we add noise, then the variance of the function in between data points becomes very large.
In the image we show this by computing the standard deviation for the sample distribution of the estimate when there is white noise with variance of $\sigma = 0.2, 0.4 \text{ or } 0.6$
In this case with fitting a 10h order polynomial you see that the error due to sample variation is relatively constant over all values for $x$. It is only for extrapolation that there is a large influence.
So this matter of overfitting does not seem to be much like Runge's phenomenon.


I have to explain the graphs a bit better. But, below is the R-code that allows to see what I did.
### create data
set.seed(1)
sigma = 0.2
### fine grained data
xs = seq(-0.1, 1.1, 0.001)
ys = sin(xs^2*6*pi)
### 21 data points
x = seq(0,1,0.05)
y = sin(x^2*6*pi)

### polynomials for modelling
M = as.data.frame(poly(x,10))
Ms = as.data.frame(predict(poly(x,10), xs))

### model of fit without noise
mod1 <- lm(y ~ ., data = M)

plot(xs,ys, type = "l", ylim = c(-2,2), lwd = 2,
     xlab = "x", ylab = "y", xlim = c(0,1))

### add confidence intervals at three levels of sigma
for (sig in 1:3) {
  pred = predict(mod1, newdata = Ms)+cor1
  ### the error of the coefficients is sigma * I because 
  ### (X^tX)^-1 is equal to the identity matrix 
  sigerr = sqrt(rowSums(Ms^2))*sigma*sig
  polygon(c(rev(xs), xs), c(rev(pred+sigerr),pred-sigerr),
          col = rgb(1,0,0,0.3), border = NA)
}

### true data
lines(xs,ys, lwd = 2)


### use this to simulate multiple fits
#for (i in 1:1000) {
#  q = y+rnorm(21,0,sigma)
#  mod <- lm(q ~ ., data = M)
#  lines(xs, predict(mod, newdata = Ms)+cor1, lty = 1, lwd = 1, col = rgb(1,0,0,0.01)) 
#}

### fit of model to 21 points
lines(xs, predict(mod1, newdata = Ms), lty = 2, lwd = 2)

## data points
points(x,y, pch = 21, col = 1, bg = 0)

title(expression(sin(6 * pi * x^2) * " with least squares fit to 21 points"))

legend(0,-1.2, c("true function with noisless data points",
                 "fit of noisless data points"), 
       lty = c(1,2), lwd = 2, pch = c(21,NA), pt.bg = c(0,0), cex = 0.8)

