Fourier transform and anti-trasform--identity missing
Different ways of showing it, depends where you start. Are you willing to accept that the Fourier transform of $\delta(t)$ is $1$ and vice versa (i.e., $\int_{-\infty}^\infty {\color{red}1} \cdot {\rm e}^{\jmath 2\pi f t} = \delta(t)$)? If so, it's easy:

$$\begin{align}
r(\tau) &amp; = \int_{-\infty}^\infty {\rm e}^{\jmath 2\pi f \tau} S(f) {\rm d}f \\
&amp; = \int_{-\infty}^\infty {\rm e}^{\jmath 2\pi f \tau} \int_{-\infty}^\infty {\rm e}^{-\jmath 2\pi f t} r(t) {\rm d}t {\rm d}f \\
&amp; = \int_{-\infty}^\infty \int_{-\infty}^\infty {\rm e}^{\jmath 2\pi f (\tau-t)} r(t) {\rm d}t {\rm d}f \\
&amp; = \int_{-\infty}^\infty r(t) \int_{-\infty}^\infty {\color{red}1} \cdot {\rm e}^{\jmath 2\pi f (\tau-t)} {\rm d}f {\rm d} t \\
&amp; = \int_{-\infty}^\infty r(t) \delta(\tau-t) {\rm d} t \\
&amp; = \int_{-\infty}^\infty r(\tau) \delta(\tau-t) {\rm d} t \\
&amp; = r(\tau) \int_{-\infty}^\infty \delta(\tau-t) {\rm d} t \\
&amp; = r(\tau)
\end{align}$$


Step 1: Insert (2) into (1). Note that the inner integration variable is a new one, different from $\tau$. I call it $t$.
Step 2: Pulling the first exp inside the integral, combining the exps.
Step 3: Changing integration order (PSD and ACF are absolutely integrable), pulling out what does not depend on the inner integration variable $f$.
Step 4: Using the fact that the inverse Fourier of a constant is a delta (think $\tau-t$ as one variable here, then it's the inverse Fourier of 1).
Step 5: Using the sifting property of the delta.
Step 6: Moving out the constant term
Step 7: Area under the delta is one.


Of course, step 4 is the critical one. If you don't buy it, you need a different, more fundamental/mathematical approach. This is more the engineering point of view I'm presenting here.

Regarding your reply with $\tau$ vs. $t$: What you say is not true. See, we're computing $r(\tau)$ via $\int_{-\infty}^\infty {\rm e}^{\jmath 2\pi f\tau} S(f) {\rm d}f$, which means it's an integral over frequency and the integration kernel depends on $\tau$. The function $S(f)$ does not depend on $\tau$ of course. Now, you are replacing $S(f)$ by the inverse Fourier transform of the autocorrelation, which is an integral over time. But it's a different time variable, which I called $t$. It must be, since it if were $\tau$ it would mean that $S(f)$ somehow depends on $\tau$. 

Another way to think about it: The variable $\tau$ is the independent one. Don't forget that all our integrals are definite ones (from $-\infty$ to $\infty$). We sometimes drop that for laziness, but I added them now to be more clear. Now, the integration variables on the right-hand side are the ones we integrate over, hence the result cannot depend on it. Imagine something like $x(\tau) = \int_{-\infty}^\infty \int_{-\infty}^\infty X(f) g(\tau) {\rm d}f {\rm d}\tau$. This does not make sense as the right-hand side integrates over $f$ and $\tau$ (the result is a number) whereas the right-hand side depends on $\tau$. This is why we need a new time variable.