Prediction intervals for simple/baseline forecasts I'm reading this, and I don't understand how the prediction intervals are calculated for the baseline forecast methods. I agree that we can estimate the SD of the sample from the known residuals
$$ \begin{equation}
  \hat{\sigma} = \sqrt{\frac{1}{T-K-M}\sum_{t=1}^T e_t^2}, \tag{5.1}
\end{equation}$$
Where $K$ is number of parameters, and $M$ is number of missing values. And I agree that intuitively they should grow with longer forecast (grow with $h$). But I'm not sure how you actually arrive to these quantities? (Assuming normally distributed and uncorrelated residuals).

Can anyone show the math behind this?

We assume that the residuals from the method are uncorrelated and homoscedastic, with mean 0 and variance $\sigma^2$. We don't need to assume normality. Let $y_1,\dots,y_T$ denote the time series observations, and let $\hat{y}_{T+h|T}$ be the estimated forecast mean (or point forecast). Then we can write
$$
y_t = \hat{y}_{t|t-1} + \varepsilon_t
$$
where $\varepsilon$ is a white noise process. Let $\hat\sigma_h^2$ be the estimated $h$-step forecast variance.
Random walk/Naive method
For a random walk, the equation above suggests that the appropriate model is
$$y_t = y_{t-1} + \varepsilon_t.$$
Therefore
\begin{align*}
y_{T+1} &= y_T + \varepsilon_T \\
y_{T+2} &= y_{T+1} + \varepsilon_{T+1}  = y_T + \varepsilon_T + \varepsilon_{T+1}\\
\vdots \\
y_{T+h} &= y_T + \sum_{i=1}^h \varepsilon_{T+i},
\end{align*}
and so
\begin{align*}
\hat{y}_{T+h|T} &= E(y_{T+h|T}) = y_T\\
\text{and}\qquad
\hat{\sigma}^2_{h} &= V(y_{T+h|T}) = h\sigma^2
\end{align*}
Seasonal naive method
Here the model is
$$y_t = y_{t-m} + \varepsilon_t,$$
where $m$ is the seasonal period. Thus
\begin{align*}
y_{T+1} &= y_{T+1-m} + \varepsilon_{T+1} \\
y_{T+2} &= y_{T+2-m} + \varepsilon_{T+2} \\
\vdots \\
y_{T+m} &= y_{T} + \varepsilon_{T+m} \\
y_{T+m+1} &= y_{T+1} + \varepsilon_{T+m+1} = y_{T+1-m} + \varepsilon_{T+1} + \varepsilon_{T+m+1}\\
\vdots\\
y_{T+2m} &= y_{T+m} + \varepsilon_{T+2m} = y_{T} + \varepsilon_{T+m} + \varepsilon_{T+2m}\\
y_{T+2m+1} &= y_{T+m+1} + \varepsilon_{T+2m+1} = y_{T+1-m} + \varepsilon_{T+1} + \varepsilon_{T+m+1} + \varepsilon_{T+2m+1}\\
\vdots\\
y_{T+h} &= y_{T+h-m(k+1)} + \sum_{i=0}^{k} \varepsilon_{T+h+m(i-k)},
\end{align*}
where $k$ is the integer part of $(h-1) /m$ (i.e., the number of complete years in the forecast period prior to time $T+h$).
Therefore
\begin{align*}
\hat{y}_{T+h|T} &= E(y_{T+h|T}) = y_{T+h-m(k+1)}\\
\text{and}\qquad
\hat{\sigma}^2_{h} &= V(y_{T+h|T}) = (k+1)\sigma^2
\end{align*}
Mean method
The model underpinning the mean method is
$$y_t = c + \varepsilon_t$$
for some constant $c$ to be estimated. The least-squares estimate of $c$ is the mean, $\hat{c} = \bar{y} = (y_1+\dots+y_T)/T.$ Thus,
\begin{align*}
y_{T+1} &= c + \varepsilon_{T+1} \\
y_{T+2} &= c + \varepsilon_{T+1} \\
\vdots \\
y_{T+h} &= c + \varepsilon_{T+h}.
\end{align*}
Therefore
\begin{align*}
\hat{y}_{T+h|T} &= E(y_{T+h|T}) = \hat{c} = \bar{y}\\
\text{and}\qquad
\hat{\sigma}^2_{h} &= V(y_{T+h|T}) = V(\hat{c}) + \sigma^2 = (T^{-1} + 1)\sigma^2
\end{align*}
Drift method
For a random walk with drift
$$y_{t} = c + y_{t-1} + \varepsilon_t.$$
Therefore,
\begin{align*}
y_{T+1} &= c + y_T + \varepsilon_T \\
y_{T+2} &= 2 + y_{T+1} + \varepsilon_{T+1}  = 2c + y_T + \varepsilon_T + \varepsilon_{T+1}\\
\vdots \\
y_{T+h} &= hc + y_T + \sum_{i=1}^h \varepsilon_{T+i}.
\end{align*}
Now the least squares estimate of $c$ is $\hat{c} = (y_{T}-y_1)/(T-1)$.
Therefore
\begin{align*}
\hat{y}_{T+h|T} &= E(y_{T+h|T}) = h\hat{c} + y_T\\
\text{and}\qquad
\hat{\sigma}^2_{h} &= V(y_{T+h|T}) = V(\hat{c}) + h\sigma^2 = \frac{h^2\sigma^2}{T-1} + h\sigma^2.
\end{align*}
I've also cross-posted this answer on my blog at https://robjhyndman.com/hyndsight/forecasting_variances.html
