How can I generate correlated timeseries made up of 0s and 1s? I want to generate series of 0s and 1s that exhibit some clustering. By this I mean that 1s and 0s should occur together. So I envisage series of 0s and 1s that will exhibit similar clustering of these elements, and not just random series of 0s and 1s.
In essence, for a single time series, I would go about that by thresholding a Markov chain with a 2x2 transition matrix, with some stochastic added on to it. Now, I'm not too certain on how to do this, but since I would like to produce several of these series, I was wondering whether there's something straightforward that I have missed.
I plan to use these series to simulate data availability (0 or 1) in a data acquisition system and to do some Monte Carlo simulations of how this affects what we can do with the data.
In order for the above simulations to be realistic, I would like to fit real observations to this model, so that the simulated time series share temporal correlation with the data. I would initially do this by calculating lag autocorrelations of both series and tweaking model parameters until I get something that resembles my observations, but unsure whether this is the best way.
Thanks!
 A: A standard method is to begin by generating an autocorrelated Gaussian process $z_i$.  (It doesn't have to be Gaussian, but such processes are easy to generate.)  Take the logistic (inverse logit) of the values, producing a series of numbers $p_i = 1/\left(1 + \exp(-z_i)\right)$ in the interval $(0,1)$.  Independently draw values from Bernoulli($p_i$) distributions to create a series of $0$ and $1$ values.  Clustering will tend to occur with positive autocorrelation.
As a bonus, this procedure allows you to perform two stages of simulation: you can fix the underlying realization of the Gaussian process and iterate the second stage of Bernoulli draws.  Or you can generate a separate realization of the Gaussian process each time.
There are probably R packages to do all this directly.  The geoRGLM package performs this simulation in two dimensions (using a Matern autocorrelation function, which includes exponential and Gaussian autocorrelations as special cases); you could simulate along a straight line (or $1$ by $n$ grid) to obtain a time series. 
