Algorithm to produce autocorrelated uniformly distributed number I would like to produce a time-series of autocorrelated probabilities (with a predefined mean level of autocorrelation). 
I've spotted this and this which I believe should give me what I'm looking for, but I'm unable to access them. 
If anyone could give me some hint on this 'sum of uniforms' method it would be appreciated. Many thanks.
 A: The first abstract might as well be the entire paper, because it perfectly describes an algorithm. :-)  The following is my interpretation of that abstract.
The idea emulates that of generating a Gaussian process with specified autocorrelation $\rho = \cos(\theta)$ (re-expressing $\rho$ as the cosine of an angle between $0$ and $\pi$). By rescaling and shifting the mean, it suffices to generate a process $X_n$ with Standard Normal marginals (that is, with zero mean and unit variance).  If we generate an iid sequence of Standard Normals $\varepsilon_n$, then the desired process is obtained recursively via
$$\eqalign{
X_0 &= \varepsilon_0 \\
X_{n+1} &= \cos(\theta)X_n + \sin(\theta)\varepsilon_{n+1}, \quad n \ge 0.
}$$
This works because


*

*$X_0$ is Standard Normal by construction,

*All the $X_n$ have Normal distributions because they are (inductively) linear combinations of independent Normals,

*$\text{Var}(X_{n+1})$ = $\cos(\theta)^2 \text{Var}(X_n) + 2\cos(\theta)\sin(\theta) \text{Cov}(X_n, \varepsilon_{n+1}) + \sin(\theta)^2 \text{Var}(\varepsilon_{n+1})$ = $\cos(\theta)^2 1 + \sin(\theta)^2 1 = 1$ by induction on $n$.

*$\text{Cov}(X_n, X_{n+1})$ = $\text{Cov}(X_n, \cos(\theta)X_n + \sin(\theta)\varepsilon_{n+1})$ = $\cos(\theta)\text{Var}(X_n, X_n) = \cos(\theta)$, as intended.
The analog of a Standard Normal distribution might be called the "standard uniform" distribution; an easy calculation establishes that it is the uniform distribution on $[-\sqrt{3}, \sqrt{3}]$.  If we try the same construction, beginning with a sequence $(\varepsilon_n)$ of iid standard uniforms, then all the variance and covariance calculations remain correct, but unfortunately, $Y_{n+1} = \cos(\theta)X_n + \sin(\theta)\varepsilon_{n+1}$ is not uniform.  We can make it uniform through the probability transform: simply apply the CDF of $Y_{n+1}$ to the value which is obtained--which will produce a uniform distribution on $[0,1]$--and then standardize that to make it standard uniform.  Although the correlation is no longer exactly $\cos(\theta)$, it turns out it will be extremely close to $\cos(\theta)$.
A simple way to compute the CDF of $Y_{n+1}$, which by definition is the CDF of a linear combination of two standard uniforms with coefficients $\cos(\theta)$ and $\sin(\theta)$, begins with the case $0 \le \theta \le \pi/4$.  Computing the integrals (which involve areas of triangles and parallelograms) gives
$$
\eqalign {
 1,\quad & z\geq \sqrt{3} (c+s) \\
 \frac{1}{12} \csc(2 \theta) \left(-3-z^2+2 \sqrt{3} z (c+s)+9 \sin(2 \theta)\right),\quad & \sqrt{3} (c-s)<z<\sqrt{3} (c+s) \\
 \frac{1}{6} \left(3+\sqrt{3} z /c\right),\quad & \sqrt{3} (-c+s)<z\leq \sqrt{3} (c-s) \\
 \frac{1}{12} \csc(2 \theta) \left(z^2+2 \sqrt{3} z (c+s)+3 (1+\sin(2 \theta))\right),\quad & -\sqrt{3} (c+s)<z\leq \sqrt{3} (-c+s) \\
 0, \quad &\text{otherwise}
}
$$
where $c=\cos(\theta)$ and $s=\sin(\theta)$.  (Computations are made a little easier by noting the distribution is symmetric.)  Let's call this function $F(z;\theta)$. Here is a plot of $F$ and its derivative (the PDF) for $\theta=1/2$:

I have broken the graph of $F$ (in blue) at the points where its definition changes.  These correspond to changes in slope of its derivative (in red, shaded).
Now exploit the symmetries available via $x\to -x$, $y \to -y$, and $(x,y)\to (y,x)$ to reduce all computations to the case $0 \le \theta \lt \pi/4$: define 
$$\psi = |\text{mod}(\text{mod}(\theta, \pi), \pi/2, -\pi/4))|$$
to be the absolute difference between $\theta$ (reduced modulo $\pi$) and the nearest multiple of $\pi/2$ and set
$$h(x,y;\theta) = 2\sqrt{3} (F(\cos(\theta)x + \sin(\theta)y; \psi) - 1/2).$$
The method of obtaining $X_{n+1}$ is to compute $h(X_n, \varepsilon_{n+1})$.

As a check, here are results obtained by simulating $10^5$ iid standard uniform $\varepsilon_n$ and using them to generate realizations of four processes with a wide scope of intended correlations at $\rho=0.9$, $0.7$, $0$, and $-0.99$, respectively.  
Their means are $0.0146058,0.00813056,0.0040066,0.000269295$.  These are close to zero and can be attributed to random variation; the large serial correlation in the first accounts for the relatively large departure of its mean from $0$.  Likewise, their variances of $1.00962,1.00654,0.99683,0.996235$ are only insignificantly different from $1$.
We visually check the distributions are uniform by plotting their histograms:

They look close to uniform, as intended.  (This really wasn't in doubt, because the construction--assuming the calculation is carried out correctly--forces them to be uniform; the main concern is whether the serial correlation coefficients are close to the intended values.)  With $10^5$ iterations we should be able to pin down the correlations to better than two decimal places.  The lag-one serial correlations are 
$$0.906885, 0.696603, 0.00768022, -0.990412,$$
comfortably close to the targets.  For further confirmation, here are plots of the first $300$ values in each realization:

