$\newcommand{\d}{\,\mathrm{d}}\newcommand{\p}{\operatorname{Pr}}\newcommand{\pdf}{\operatorname{pdf}}\newcommand{\cdf}{\operatorname{cdf}}\newcommand{\n}{\mathcal{N}}$Suppose we have a sequence of $n$ standard normal variables $\{X_k\}_{k=1}^n$ which are i.i.d: I want to show that:
$$\sum_{k=1}^n(X_k-\bar{X})^2\sim\chi^2_{n-1}$$
First, I use these formulae, which can be derived by differentiating the $\cdf$ for independent variables:
$$\pdf_{X^2}(x)=\begin{cases}\frac{1}{2\sqrt{x}}\pdf_X(\sqrt{x})-\frac{1}{2\sqrt{x}}\pdf_X(-\sqrt{x})&x\gt 0\\0&x\le 0\end{cases}\\\pdf_{X+Y}(x)=(\pdf_X\ast\pdf_Y)(x)\\\pdf_{\alpha X}(x)=\frac{1}{|\alpha|}\pdf_X(x/\alpha)$$
Now I expand the sum, which is a random variable I'll call $S$:
$$S=\color{red}{\sum_{k=1}^nX^2_k}\color{green}{-\frac{1}{n}\left(\color{blue}{\sum_{k=1}^nX_k}\right)^2}$$
The red, green and blue random variables I'll call $A,B,C$.
I show that $A$ is distributed according to $\chi^2_n$:
First let us check $n=1$, i.e. the distribution of a squared standard normal variable $X$. According to the formula, when $x\gt0$: $$\pdf_X^2(x)=\frac{1}{2\sqrt{x}}\pdf_X(x)-\frac{1}{2\sqrt{x}}\pdf_X(-\sqrt{x})=\frac{x^{-1/2}}{\sqrt{2\pi}}\exp\left(-\frac{1}{2}x\right)$$Which is the $\pdf$ of $\chi^2_1$ as expected. Now inductively assume that $\sum_{k=1}^{n}X_k\sim\chi^2_n$, and consider $\sum_{k=1}^{n+1}X_k^2$, where $X_k$ are again all standard normal variables. Using the addition formula, the $\pdf$ of the new sum will be the convolution between two $\chi^2$ distributions, as the additional $X_{n+1}^2$ was shown to follow $\chi^2_1$: $$\pdf_{\sum_{k=1}^{n+1}X^2_k}(x)=\int_{-\infty}^\infty\pdf_{\chi^2_n}(t)\pdf_{\chi^2_1}(x-t)\d t$$But since the $\chi^2$ $\pdf$ is zero for negative values, the integral becomes $0$ for negative $x$, and for positive $x$ it becomes: $$\int_0^x\left(\frac{t^{n/2-1}\exp\left(-\frac{1}{2}t\right)}{2^{n/2}\Gamma(n/2)}\right)\left(\frac{(x-t)^{-1/2}\exp\left(-\frac{1}{2}(x-t)\right)}{\sqrt{2\pi}}\right)\d t$$Which can be simplified to: $$\frac{\exp\left(-\frac{1}{2}x\right)}{2^{(n+1)/2}\Gamma(n/2)\Gamma(1/2)}\int_0^xt^{n/2-1}(x-t)^{-1/2}\d t$$Under $u=t/x$, we find: $$\frac{\exp\left(-\frac{1}{2}x\right)}{2^{(n+1)/2)}\Gamma(n/2)\Gamma(1/2)}\int_0^1(ux)^{n/2-1}(x-ux)^{-1/2}(x\d u)=\frac{x^{(n+1)/2-1}\exp\left(-\frac{1}{2}x\right)}{2^{(n+1)/2)}\Gamma(n/2)\Gamma(1/2)}\mathcal{B}\left(\frac{n}{2},\frac{1}{2}\right)$$Finally giving: $$\pdf_{\sum_{k=1}^{n+1}X_k^2}(x)=\frac{x^{(n+1)/2-1}\exp\left(-\frac{1}{2}x\right)}{2^{(n+1)/2}\Gamma(n/2)\Gamma(1/2)}\cdot\frac{\Gamma(n/2)\Gamma(1/2)}{\Gamma((n+1)/2)}=\pdf_{\chi^2_{n+1}}(x)$$Inductively therefore we are done, and the variable $A$ follows a $\chi^2_n$ distribution.
Now examine the blue variable, $C$:
I will show inductively that: $$\pdf_{\sum_{k=1}^nX_k}(x)=\frac{\exp\left(-\frac{1}{2n}x^2\right)}{\sqrt{2\pi\cdot n}}$$This is clearly true of $n=1$. Assume its truth for some $n$ and consider the case $n+1$; adding by convolution gives: $$\begin{align}\pdf_{\sum_{k=1}^{n+1}X_k}(x)&=\frac{\exp\left(-\frac{1}{2}x^2\right)}{\sqrt{2\pi\cdot n}\sqrt{2\pi}}\int_{-\infty}^\infty\exp\left(-\frac{1}{2}(t^2-2t\cdot x)-\frac{1}{2n}t^2\right)\d t\\&=\frac{\exp\left(-\frac{1}{2}x^2\right)}{\sqrt{2\pi\cdot n}\sqrt{2\pi}}\int_{-\infty}^\infty\exp\left(-\left(\frac{1}{2}+\frac{1}{2n}\right)\left(t^2-\left(\frac{1}{2}+\frac{1}{2n}\right)^{-1}t\cdot x\right)\right)\d t\\&=\frac{\exp\left(-\frac{1}{2}x^2\right)}{\sqrt{2\pi\cdot n}\sqrt{2\pi}}\int_{-\infty}^\infty\exp\left(-\frac{n+1}{2n}\left(t^2-\frac{2n}{n+1}t\cdot x\right)\right)\d t\\&=\frac{\exp\left(-\frac{1}{2}x^2\right)}{\sqrt{2\pi\cdot n}\sqrt{2\pi}}\int_{-\infty}^\infty\exp\left(-\frac{n+1}{2n}\left(t-\frac{n}{n+1}\cdot x\right)^2+\frac{n+1}{2n}\left(\frac{n}{n+1}\right)^2x^2\right)\d t\\&=\frac{1}{\sqrt{2\pi\cdot n}\sqrt{2\pi}}\exp\left(\left(\frac{n}{2}\frac{1}{n+1}-\frac{1}{2}\right)x^2\right)\int_{-\infty}^{\infty}\exp\left(-\frac{n+1}{2n}t^2\right)\d t\\&=\frac{1}{\sqrt{2\pi\cdot n}\sqrt{2\pi}}\exp\left(-\frac{1}{2(n+1)}x^2\right)\sqrt{\frac{2\pi\cdot n}{n+1}}\\&=\frac{1}{\sqrt{2\pi(n+1)}}\exp\left(-\frac{1}{2(n+1)}x^2\right)\end{align}$$Inductively, we are done.
Now let's look at the green; $C^2$, which is then scaled by $-\frac{1}{n}$:
Using the evenness of $C$, $$\pdf_{C^2}(x)=\frac{1}{\sqrt{x}}\pdf_C(\sqrt{x})=\frac{1}{\sqrt{2\pi nx}}\exp\left(-\frac{1}{2n}x\right)$$Scaling by $-\frac{1}{n}$ yields the green variable, $B$, to have $\pdf$, now only valid for negative $x$: $$\pdf_{B}(x)=\frac{1}{\sqrt{2\pi}}\cdot|x|^{-1/2}\exp\left(\frac{1}{2}x\right)\sim-\chi^2_1$$And the $\pdf$ is now $0$ for positive $x$.
Finally, $S$ has distribution equivalent to $A+B$, and so we convolve:
$$\pdf_S(x)=\int_{-\infty}^{\infty}\pdf_A(t)\pdf_B(x-t)\d t$$
But numerical evaluation suggests that this last integral is not the desired $\pdf$ for $\chi^2_{n-1}$. What's my mistake?
Note: I realise that $B$ is simply the negative of a $\chi^2_1$ variable, and I could try to argue that since $\chi^2_{n-1}+\chi^2_1\sim\chi^2_n$, then $\chi^2_n-\chi^2_1\sim\chi^2_{n-1}$. However, I believe that this is not true, since the subtraction of two independent normal variables yields the same as the addition, and the variables are independent, right? Maybe the solution is that since $A,B$ both depend on $\{X_k\}_{k=1}^n$, they are not independent so the convolution formula doesn't hold, and indeed I can make the seemingly naive subtraction of $\chi^2_n-\chi^2_1\sim\chi^2_{n-1}$
Best Answer
We first note that, if $\mathbf{1} = (1, \ldots, 1)$ and $\mathrm{X} = (X_1,\ldots,X_n)$ are regarded as column vectors, then $S$ is a quadratic form in $X$:
$$ S = \sum_{k=1}^{n} (X_k - \overline{X})^2 = \sum_{k=1}^{n} X_k^2 - \frac{1}{n} (X^{\mathsf{T}} \mathbf{1})(\mathbf{1}^{\mathsf{T}} X) = X^{\mathsf{T}} P X, $$
where $P = I_n - \frac{1}{n}\mathbf{1}\mathbf{1}^{\mathsf{T}}$ is the orthogonal projection onto the orthogonal complement $\operatorname{span}(\mathbf{1})^{\perp}$. So, the moment generating function of $S$ is given by
\begin{align*} \mathbf{E}[e^{\lambda S}] = \int_{\mathbb{R}^n} e^{\lambda x^{\mathsf{T}}Px} \frac{e^{-|x|^2/2}}{(2\pi)^{n/2}} \, \mathrm{d}x = \frac{1}{(2\pi)^{n/2}} \int_{\mathbb{R}^n} e^{-x^{\mathsf{T}}\left( \frac{1}{2}I_n - \lambda P \right)x} \, \mathrm{d}x. \end{align*}
Now if $\Sigma$ is an $n\times n$ positive-definite matrix, then it is well-known that
$$ \int_{\mathbb{R}^n} e^{-x^{\mathsf{T}}\Sigma x} \, \mathrm{d}x = \frac{\pi^{n/2}}{\sqrt{\det \Sigma}}. $$
Since the eigenvalues of $\frac{1}{2}I_n - \lambda P$ are $\frac{1}{2}, \frac{1}{2}-\lambda, \ldots, \frac{1}{2}-\lambda$, it follows that $\frac{1}{2}I_n - \lambda P$ is positive-definite for $\lambda < \frac{1}{2}$ and
\begin{align*} \mathbf{E}[e^{\lambda S}] &= \frac{1}{2^{n/2}\sqrt{\det\left( \frac{1}{2}I_n - \lambda P \right)}} = \frac{1}{\sqrt{\det\left( I_n - 2 \lambda P \right)}} = \frac{1}{(1-2\lambda)^{(n-1)/2}}. \end{align*}
This is precisely the m.g.f. of $\chi^2_{n-1}$, and therefore $S \sim \chi^2_{n-1}$.