The other answer doesn't really give the full story. It's the transformation theorem. I just had an answer involving it here.
Call
$$X(t_1) = Y_1,~X(t_2) = Y_2 + Y_1, ...,~X(t_n) = Y_n + \cdots + Y_1.$$
The inverse transformations are $$Y_1 = X(t_1),~Y_2 = X(t_2)-X(t_1), ..., ~Y_n = X(t_n) - X(t_{n-1}).$$
We know the joint density of all the $Y$s is the product of mean zero gaussians with variances $\sigma^2 t_1$, $\sigma^2 (t_2-t_1),\ldots, \sigma^2(t_n-t_{n-1})$. We're using the independent and stationary increments properties here. We will call this product of Gaussians $f_{Y_1,\ldots,Y_n}(y_1,\ldots,y_n) = f_{Y_1}(y_1)\cdots f_{Y_n}(y_n)$. Looking at the matrix of partial derivatives, its determinant is $1$.
So this is why you just plug stuff in. The density of $X(t_1),\ldots, X(t_n)$ is
\begin{align*}
g_{X(t_1),\ldots, X(t_n)}(x_1,\ldots,x_n)&=f_{Y_1,\ldots,Y_n}(y_1[x(t_1),\ldots,x(t_n)],\ldots,y_n[x(t_1),\ldots,x(t_n)])|1|\\
&= f_{Y_1,\ldots,Y_n}(x(t_1),x(t_2)-x(t_1),\ldots,x(t_n)-x(t_{n-1})) \\
&= f_{Y_1}(x(t_1))f_{Y_2}(x(t_2)-x(t_1))\cdots f_{Y_n}(x(t_n)-x(t_{n-1})).
\end{align*}
Best Answer
If you don't worry about technicality too much, there's a nice proof on Wikipedia. The gist of it is as follows: we know that white noise $\xi(t)$ has a flat power spectrum (just simulate some to see this!) and also that the Wiener process is the integral of white noise, so that we must have $$ S_{xx}^{(\xi)}(\omega) = \left| \mathcal{F}\left\{ \frac{d\mathcal{W}}{dt} \right\} \right|^2 \propto 1. $$ Recalling that $\mathcal{F}\left\{ \frac{d}{dt}f(t) \right\} \propto \omega\hat{f}(\omega)$, we see that we must have $\hat{f}(\omega) \propto \frac{1}{\omega}$. Then $S_{xx}^{(\mathcal{W})}(\omega) = \left| \hat{f}(\omega) \right|^2 \propto \frac{1}{\omega^2},$ completing the proof. Now, if you are concerned about some technicalities of convergence, there is no problem either. Here is what you can do. Redefine the Fourier transform as $$ \mathcal{F}\left\{f(t)\right\} = \frac{1}{\sqrt{T}}\int_0^T f(t) \exp(-i \omega t)\ dt $$ and note that this mean exists for all finite $T$. Then define the power spectral density as $$ S_{xx}(\omega) = \lim_{T \rightarrow \infty}\left \langle \left| \mathcal{F}\{x(t)\}\right|^2 \right\rangle, $$ where the mean is over realizations of the process $x(t)$.