Why do the widths of confidence & prediction intervals change across a regression line - shouldn't it be the same with i.i.d? I thought that we assume i.i.d and therefore identically distributed (constant variance) over the line? But aren't the below saying that technically the distribution changes from one X value to another?
Example images of confidence intervals:


I thought that we assume i.i.d and therefore identically distributed (constant variance) over the line? But aren't the above saying that technically the distribution changes from one X value to another?
Edit
I think what I'm saying may be better described by the below image.
If we had enough data and could theoretically make a slice at each value of x_1, so as to obtain a theoretical conditional probability distribution of x_2 | x_1, wouldn't we expect this distribution to be the same around the best fit line for all x_1, for it to satisfy i.i.d? In the case of say normally distributed prediction intervals that widen as we go further from the centre of the data, it seems to me that we're suggesting the variance is not constant and therefore the distributions are not identical for all x.
Sorry if this makes no sense.


I think i can give a more mathematical answer then the two already provided.
You are fitting a linear model $Y = X\beta + \epsilon$. $X$ being a design matrix  $\epsilon \sim N(0, \sigma_\epsilon)$.
The estimators of your model parameters $\hat\beta$ have an estimated covariance matrix $\hat{\Sigma}_\beta$, from which the standard error of  $ \hat E(Y|X=x) = \hat y(x) = x^T\hat\beta$(with $X$ now being a general stand in for the independent variables) can be calculated with
$$
SE(\hat y(x)) = \sqrt{x^T\hat{\Sigma}_\beta x}
$$
This will almost never be a constant, so the width of your confidence interval will vary.
Now for the predictive interval you assumed that $(Y|X=x) \sim N(\hat y(x), \hat\sigma_\epsilon)$, but you should have included the uncertainty about $\hat y(x)$, so $(Y|X=x) \sim \text{Student_t}\left(\mu = \hat y(x), \sigma =\sqrt{\hat\sigma_\epsilon^2 + SE(\hat y(x))^2}, \nu = \text{model residual df} \right)$
Here is some R code to illustrate the point:
n <- 50
x <- rnorm(n)
y <- 0.5*x + rnorm(n, sd = 1) 

m <- lm(y ~x)
x_help <- seq(-4, 4, length.out = 1000)

sum_m <- summary(m)
res.sig <-sum_m$sigma
df_conf <-  as.data.frame(predict(m, newdata = data.frame(x = x_help), interval = "conf", se.fit = T))
df_pred <- as.data.frame(predict(m, newdata = data.frame(x = x_help), interval = "predict"))
plot(x_help, df_pred$upr - df_pred$lwr, 
     ylim = c(2*qnorm(0.975)*res.sig, max( df_pred$upr - df_pred$lwr)),
     ylab = "Width of 95%-prediction intervall")
abline(v = mean(x))
# naive straight 
abline(h = 2*qnorm(0.975)*res.sig)
# adding variances, but normal distribution is wrong
lines(x_help, 2*qnorm(0.975)*sqrt(res.sig^2 + df_conf$se.fit^2))
# using t-distribution
lines(x_help, 2*qt(0.975, df = n - 2)*sqrt(res.sig^2 + df_conf$se.fit^2), col = 2) 

and here is the resulting figure: 
