Como gerar dados binários aleatórios de séries temporais correlacionados automaticamente?

15

Como posso gerar séries temporais binárias, de modo que:

  1. A probabilidade média de observar 1 é especificada (digamos 5%);
  2. Probabilidade condicional de observar 1 no tempo t dado o valor em (digamos 30% se o valor de for 1)?t1t1
user333
fonte

Respostas:

17

Use uma cadeia de Markov de dois estados.

Se os estados são chamados 0 e 1, a cadeia pode ser representada por uma matriz 2x2 fornecendo as probabilidades de transição entre estados, onde P i j é a probabilidade de passar do estado i para o estado j . Nesta matriz, cada linha deve somar 1.0.PPijij

Na afirmação 2, temos , e a conservação simples diz P 10 = 0,7 .P11=0.3P10=0.7

Na afirmação 1, você deseja que a probabilidade de longo prazo (também chamada de equilíbrio ou estado estacionário) seja . Isso indica que P 1 = 0,05 = 0,3 P 1 + P 01 ( 1 - P 1 ) A solução fornece P 01 = 0,0368421 e uma matriz de transição P = ( 0,963158 0,0368421 0,7 0,3 )P1=0.05

P1=0.05=0.3P1+P01(1P1)
P01=0.0368421
P=(0.9631580.03684210.70.3)

(Você pode verificar se a sua matriz de transição está correta, elevando-a a uma potência alta - nesse caso, 14 faz o trabalho - cada linha do resultado fornece as mesmas probabilidades de estado estacionário)

Agora, no seu programa de números aleatórios, comece escolhendo aleatoriamente o estado 0 ou 1; isso seleciona qual linha de você está usando. Em seguida, use um número aleatório uniforme para determinar o próximo estado. Cuspa esse número, enxágue, repita conforme necessário.P

Mike Anderson
fonte
Solução interessante! Você talvez tenha algum código de exemplo no R? Antone mais?
User333
@ Mike Você pode registrar sua conta? Você é um usuário bastante ativo e precisamos mesclá-lo manualmente uma e outra vez. O processo é bastante fácil; basta visitar stats.stackexchange.com/login
Obrigado. Como posso estimar a cadeia de Markov (matriz de transição) dados os dados? Existe uma função R para fazer isso?
user333
6

Fiz uma rachadura na codificação da resposta de Mike Anderson em R. Não consegui descobrir como fazê-lo usando sapply, então usei um loop. Alterei os probs ligeiramente para obter um resultado mais interessante e usei 'A' e 'B' para representar os estados. Diz-me o que pensas.

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: Em resposta ao comentário de Paulo, aqui está uma formulação mais 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))

Eu escrevi o código original quando estava aprendendo R, então me dê um pouco de folga. ;-)

Veja como você estimaria a matriz de transição, considerando a série:

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

A ordem é trocada em relação à minha matriz de transição original, mas obtém as probabilidades corretas.

Zach
fonte
Ótimo! Eu vou assim que pissible ... Parece bom o suficiente ....
user333
É possível fazer o inverso? Dada a série estimar a matriz?
user333
Pr(Xt=Eu|Xt-1=j)
+1, mas também tenho alguns comentários: um forloop seria um pouco mais limpo aqui, você sabe o tamanho Series, então use for(i in 2:length(Series)). Isso elimina a necessidade de i = i + 1. Além disso, por que primeiro exemplo Ae depois converter para 0,1? Você pode provar diretamente 0's e 1' s.
Paul Hiemstra
2
Mais geralmente você poderia então envolvê-la em uma nova função 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);para permitir arbitrária, pré-especificado lag 1 autocorrelação
Tom Wenseleers
1

Aqui está uma resposta baseada no markovchainpacote que pode ser generalizado para estruturas de dependência mais complexas.

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)

Isso lhe dá:

enter image description here

tchakravarty
fonte
0

Perdi o controle do artigo em que essa abordagem foi descrita, mas aqui vai.

Decomponha a matriz de transição em

T=(1-pt)[10 00 01]+pt[p0 0p0 0(1-p0 0)(1-p0 0)]=(1-pt)Eu+ptE

que, intuitivamente, corresponde à ideia de que existe alguma probabilidade 1-pt que o sistema permanece no mesmo estado, e uma probabilidade pt que o estado é randomizado, onde randomizado significa fazer um sorteio independente da distribuição de equilíbrio para o próximo estado (p0 0 é a probabilidade de equilíbrio de estar no primeiro estado).

Observe que, a partir dos dados que você especificou, é necessário solucionar pt do especificado T11 através da T11=(1-pt)+pt(1-p0 0).

Uma das características úteis dessa decomposição é que ela generaliza diretamente para a classe de modelos de Markov correlacionados em problemas dimensionais mais altos.

Dave
fonte
Se alguém viu o documento que desenvolve essa representação, informe-me.
Dave