I would like to understand pink (or $1/f$) noise better. However, clearly written resources are difficult to come by, and are usually concerned with its Fourier spectrum, or qualitative descriptions of it as a time-dependent process.
I hope that someone can describe how a pink noise process $P(t)$ depends on its values at earlier times $t' < t$, in the form of a probability density function on $P(t) - P(t')$, in the same way that we can do for Brown noise (i.e. Brownian motion).
Preliminaries
As a 'process', white noise can be modelled as a collection of i.i.d. random variables $W(t)$, where for each $t$ the distribution of $W(t)$ is uniform over some interval $[-a,a]$. Brown or red noise is the time-integral of white noise, in the sense that $$ B(t) = \int_0^t W(s) \,\mathrm d s \;.$$ Note that $B(t)$ is not independent of $B(t')$ for $t' < t$, in that $B_t - B_{t'}$ is normally distributed with mean $0$ and variance proportional to $(t - t')$. Thus we can talk about the expected properties of the trajectory at discrete time intervals in a straight-forward way.
If we take the Fourier transform of a white noise process over a time interval, we obtain a function whose expected modulus-square is constant over all frequencies, albeit with increasing variance for higher frequencies. If we do the same for Brown noise, we obtain a function whose expected modulus-square varies as $1/f^2$ with respect to frequency, again with increasing variance for higher frequencies.
'Defining' pink noise
Intermediate to these is pink noise, which is a process whose Fourier transform has a modulus-square proportional to $1/f$ (so that it is often called $1/f$ noise). For a discrete-time process, where we define $W(t)$ on the integers $t \in [0,T]$ for some $T$, we may mock up a pink noise process by setting $$ \begin{aligned} \hat P(0) &= \hat W(0) \\ \hat P(k) &= \hat W(k) \big/ \sqrt{\lvert k \rvert}, \quad k \in [-T{+}1,T] \setminus \{ 0 \} \end{aligned} $$ and then defining $P$ to be the inverse Fourier transform of $\hat P$: this is what the clearest sources that I can find suggest.
A small amount of numerical experimentation on my part suggests that the resulting function $P(t)$ has a much smaller variance than the original white-noise process $W(t)$. The time-integral $C(t)$ (i.e., a cumulative sum) of $P(t)$ describes a trajectory which is less volatile than the Brown noise process $B(t)$ obtained from the same white-noise function, but otherwise the two time-integrated noise processes follow a similar trajectory. Both of these observations are not very surprising, given that $P$ arises by decreasing the importance of the high-frequency components of $W$.
Question.
There are some qualitative descriptions of pink noise, in terms of how $P(t)$ depends on $P(t')$ for $t' < t$, such as that the autocorrelation of $P(t)$ with $P(t')$ decreases exponentially towards some non-zero value and stays there. However, I would like a quantitative result similar to what we have for Brown noise — specifically: is there a precise description for the probability density function of $P(t) - P(t')$?
(Failing this, is there a better description for how to simulate a pink noise process than the one that I've outlined?)