How to generate random auto correlated binary time series data? How can I generate binary time series such that:


*

*Average probability of observing 1 is specified (say 5%);

*Conditional probability of observing 1 at time $t$ given the value at $t-1$ ( say 30% if $t-1$ value was 1)?

 A: Use a two-state Markov chain.  
If the states are called 0 and 1, then the chain can be represented by a 2x2 matrix $P$ giving the transition probabilities between states, where $P_{ij}$ is the probability of moving from state $i$ to state $j$.  In this matrix, each row should sum to 1.0.
From statement 2, we have $P_{11} = 0.3$, and simple conservation then says $P_{10} = 0.7$.
From statement 1, you want the long-term probability (also called equilibrium or steady-state) to be $P_1 = 0.05$.  This says $$P_1 = 0.05 = 0.3 P_1 + P_{01}(1-P_1)$$  Solving gives $$P_{01} = 0.0368421$$ and a transition matrix $$P = \left(
\begin{array}{cc}
 0.963158 & 0.0368421 \\
 0.7 & 0.3
\end{array}
\right)$$
(You can check your transtion matrix for correctness by raising it to a high power--in this case 14 does the job--each row of the result gives the identical steady state probabilities)
Now in your random number program, start by randomly choosing state 0 or 1; this selects which row of $P$ you're using.  Then use a uniform random number to determine the next state.  Spit out that number, rinse, repeat as necessary.
A: Here is an answer based on the markovchain package that can be generalized to more complex dependence structures.
library(markovchain)
library(dplyr)

# define the states
states_excitation = c("steady", "excited")

# transition probability matrix
tpm_excitation = matrix(
  data = c(0.2, 0.8, 0.2, 0.8), 
  byrow = TRUE, 
  nrow = 2,
  dimnames = list(states_excitation, states_excitation)
)

# markovchain object
mc_excitation = new(
  "markovchain",
  states = states_excitation,
  transitionMatrix = tpm_excitation,
  name = "Excitation Transition Model"
)

# simulate
df_excitation = data_frame(
  datetime = seq.POSIXt(as.POSIXct("01-01-2016 00:00:00", 
                                   format = "%d-%m-%Y %H:%M:%S", 
                                   tz = "UTC"), 
                        as.POSIXct("01-01-2016 23:59:00", 
                                   format = "%d-%m-%Y %H:%M:%S", 
                                   tz = "UTC"), by = "min"),
  excitation = rmarkovchain(n = 1440, mc_excitation))

# plot
df_excitation %>% 
  ggplot(aes(x = datetime, y = as.numeric(factor(excitation)))) + 
  geom_step(stat = "identity") + 
  theme_bw() + 
  scale_y_discrete(name = "State", breaks = c(1, 2), 
                   labels = states_excitation)

This gives you: 

A: I've lost track of the paper where this approach was described, but here goes.
Decompose the transition matrix into
$$
\begin{aligned}
T &= (1-p_t) \left[ \begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix} \right] + p_t \left[ \begin{matrix} p_0 & p_0 \\ (1-p_0) & (1-p_0) \end{matrix} \right] \\
&= (1-p_t) I + p_t E
\end{aligned}
$$
which, intuitively, corresponds to the idea that there is some probability $1-p_t$ that the system stays in the same state, and a probability $p_t$ that the state gets randomized, where randomized means making an independent draw from the equilibrium distribution for the next state ($p_0$ is the equilibrium probability for being in the first state).
Note that from the data you've specified you need to solve for $p_t$ from the specified $T_{11}$ via $T_{11} = (1-p_t)+p_t(1-p_0)$.
One of the useful features of this decomposition is that it pretty straightforwardly generalizes to class of correlated Markov models in higher dimensional problems.
