# How to generate random auto correlated binary time series data?

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)?

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.

• Interesting solution! Do you maybe have some sample code in R? Antone else? Aug 12, 2011 at 12:07
• @Mike Can you please register your account? You are quite an active user and we have to merge it manually over and over again. The process is quite easy; just visit stats.stackexchange.com/login
– user88
Aug 12, 2011 at 14:13
• Thanks. How can I estimate Markov chain (transition matrix) given the data? Is there an R function for doing that? Aug 12, 2011 at 21:37

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.

• Great! I'll as soon as pissible... Looks good enough.... Aug 12, 2011 at 20:16
• Is it possible to do the inverse? Given the series estimate the matrix? Aug 12, 2011 at 21:10
• "Estimate" is the key word here. The entries in the Markov chain are nothing but the conditional probabilities $Pr(X_t=i|X_{t-1}=j)$. That means you can just scan the series and find the proportions of zeroes followed by zeroes/ones and ones followed by zeroes/ones to estimate each row in the matrix. Going row by row, the usual confidence intervals obtain. Nice implementation. Aug 15, 2011 at 9:47
• +1, but I also have some comments: A for loop would be a bit cleaner here, you know the length of Series, so just use for(i in 2:length(Series)). This eliminates the need for i = i + 1. Also, why first sample A, and then convert to 0,1? You could directly sample 0's and 1's. Apr 23, 2013 at 9:15
• More generally you could then wrap it in a new function createAutocorBinSeries = function(n=100,mean=0.5,corr=0) { p01=corr*(1-mean)/mean createSeries(n,matrix(c(1-p01,p01,corr,1-corr),nrow=2,byrow=T)) };createAutocorBinSeries(n=100,mean=0.5,corr=0.9);createAutocorBinSeries(n=100,mean=0.5,corr=0.1); to allow for arbitrary, pre-specified lag 1 autocorrelation Jan 29, 2015 at 9:04

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

# 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:

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.

• If anyone has seen the paper that develops this representation, please let me know.
– Dave
Nov 21, 2018 at 15:57