How to generate random auto correlated binary time series data?

Question

How can I generate binary time series such that:

1. Average probability of observing 1 is specified (say 5%);
2. Conditional probability of observing 1 at time $t$ given the value at $t-1$ ( say 30% if $t-1$ value was 1)?

Answer 1

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.

Answer 2

I took a crack at coding @Mike Anderson answer in R. I couldn't figure out how to do it using sapply, so I used a loop. I changed the probs slightly to get a more interesting result, and I used 'A' and 'B' to represent the states. Let me know what you think.

set.seed(1234)
TransitionMatrix <- data.frame(A=c(0.9,0.7),B=c(0.1,0.3),row.names=c('A','B'))
Series <- c('A',rep(NA,99))
i <- 2
while (i <= length(Series)) {
    Series[i] <- ifelse(TransitionMatrix[Series[i-1],'A']>=runif(1),'A','B')
    i <- i+1
}
Series <- ifelse(Series=='A',1,0)
> Series
  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1
 [38] 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [75] 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1

/edit: In response to Paul's comment, here's a more elegant formulation

set.seed(1234)

createSeries <- function(n, TransitionMatrix){
  stopifnot(is.matrix(TransitionMatrix))
  stopifnot(n>0)

  Series <- c(1,rep(NA,n-1))
  random <- runif(n-1)
  for (i in 2:length(Series)){
    Series[i] <- TransitionMatrix[Series[i-1]+1,1] >= random[i-1]
  }

  return(Series)
}

createSeries(100, matrix(c(0.9,0.7,0.1,0.3), ncol=2))

I wrote the original code when I was just learning R, so cut me a little slack. ;-)

Here's how you would estimate the transition matrix, given the series:

Series <- createSeries(100000, matrix(c(0.9,0.7,0.1,0.3), ncol=2))
estimateTransMatrix <- function(Series){
  require(quantmod)
  out <- table(Lag(Series), Series)
  return(out/rowSums(out))
}
estimateTransMatrix(Series)

   Series
            0         1
  0 0.1005085 0.8994915
  1 0.2994029 0.7005971

The order is swapped vs my original transition matrix, but it gets the right probabilities.