15
$\begingroup$

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)?
$\endgroup$

4 Answers 4

17
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ Interesting solution! Do you maybe have some sample code in R? Antone else? $\endgroup$
    – user333
    Aug 12, 2011 at 12:07
  • $\begingroup$ @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 $\endgroup$
    – user88
    Aug 12, 2011 at 14:13
  • $\begingroup$ Thanks. How can I estimate Markov chain (transition matrix) given the data? Is there an R function for doing that? $\endgroup$
    – user333
    Aug 12, 2011 at 21:37
6
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ Great! I'll as soon as pissible... Looks good enough.... $\endgroup$
    – user333
    Aug 12, 2011 at 20:16
  • $\begingroup$ Is it possible to do the inverse? Given the series estimate the matrix? $\endgroup$
    – user333
    Aug 12, 2011 at 21:10
  • $\begingroup$ "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. $\endgroup$ Aug 15, 2011 at 9:47
  • $\begingroup$ +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. $\endgroup$ Apr 23, 2013 at 9:15
  • 2
    $\begingroup$ 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 $\endgroup$ Jan 29, 2015 at 9:04
1
$\begingroup$

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:

enter image description here

$\endgroup$
0
$\begingroup$

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.