¿Cómo generar datos de series de tiempo binarias correlacionadas automáticamente al azar?

¿Cómo puedo generar series de tiempo binarias tales que:

1. Se especifica la probabilidad promedio de observar 1 (digamos 5%);
2. ¿Probabilidad condicional de observar 1 en el tiempo $t$ dado el valor en (digamos 30% si el valor de era 1)?$t-1$$t-1$

Use una cadena de Markov de dos estados.

Si los estados se llaman 0 y 1, entonces la cadena se puede representar mediante una matriz 2x2 que $P$proporciona las probabilidades de transición entre estados, donde $P_{ij}$ es la probabilidad de pasar del estado $i$ al estado $j$ . En esta matriz, cada fila debe sumar 1.0.

De la declaración 2, tenemos $P_{11} = 0.3$ , y la conservación simple luego dice $P_{10} = 0.7$ .

De la declaración 1, desea que la probabilidad a largo plazo (también llamada equilibrio o estado estacionario) sea . Esto dice P 1 = 0.05 = 0.3 P 1 + P 01 ( 1 - P 1 ) Resolver da P 01 = 0.0368421 y una matriz de transición P = ( 0.963158 0.0368421 0.7 0.3$P_1 = 0.05$

(Puede verificar la corrección de su matriz de transición elevándola a una alta potencia, en este caso 14 hace el trabajo, cada fila del resultado da las mismas probabilidades de estado estacionario)

Ahora en su programa de números aleatorios, comience eligiendo aleatoriamente el estado 0 o 1; esto selecciona qué fila de estás usando. Luego use un número aleatorio uniforme para determinar el siguiente estado. Escupir ese número, enjuagar, repetir según sea necesario.$P$

Intenté codificar la respuesta de @Mike Anderson en R. No pude entender cómo hacerlo usando sapply, así que usé un bucle. Cambié los problemas ligeramente para obtener un resultado más interesante, y usé 'A' y 'B' para representar los estados. Déjame saber lo que piensas.

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: En respuesta al comentario de Paul, aquí hay una formulación más elegante

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

Escribí el código original cuando estaba aprendiendo R, así que no me preocupe. ;-)

Aquí le mostramos cómo calcularía la matriz de transición, dada la serie:

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

El orden se intercambia con mi matriz de transición original, pero obtiene las probabilidades correctas.

Aquí hay una respuesta basada en el markovchainpaquete que puede generalizarse a estructuras de dependencia más complejas.

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)

Esto te da:

He perdido la noción del documento donde se describió este enfoque, pero aquí va.

Descomponer la matriz de transición en

$1-p_t$$p_t$$p_0$

$p_t$$T_{11}$$T_{11} = (1-p_t)+p_t(1-p_0)$

Una de las características útiles de esta descomposición es que se generaliza bastante directamente a la clase de modelos de Markov correlacionados en problemas dimensionales superiores.