Wie werden zufällige automatisch korrelierte binäre Zeitreihendaten generiert?

15

Wie kann ich binäre Zeitreihen erzeugen, so dass:

  1. Die durchschnittliche Wahrscheinlichkeit, 1 zu beobachten, ist angegeben (sagen wir 5%);
  2. Bedingte Wahrscheinlichkeit, 1 zum Zeitpunkt wenn der Wert bei (sagen wir 30%, wenn der Wert von war)?t - 1 t - 1tt-1t-1
user333
quelle

Antworten:

17

Verwenden Sie eine Markov-Kette mit zwei Zuständen.

Wenn die Zustände 0 und 1 heißen, kann die Kette durch eine 2 × 2-Matrix , die die Übergangswahrscheinlichkeiten zwischen Zuständen angibt, wobei die Wahrscheinlichkeit ist, von Zustand zu Zustand . In dieser Matrix sollte jede Zeile 1,0 ergeben.P i j i jPPichjichj

Aus Aussage 2 ergibt sich , und bei einfacher Erhaltung ergibt sich dann .P 10 = 0,7P11=0,3P10=0,7

Ab Aussage 1 soll die Langzeitwahrscheinlichkeit (auch Gleichgewicht oder Steady-State genannt) . Dies besagt Das Lösen ergibt und eine ÜbergangsmatrixP 1 = 0,05 = 0,3 P 1 + P 01 ( 1 - P 1 ) P 01 = 0,0368421 P = ( 0,963158 0,0368421 0,7 0,3 )P1=0,05

P1=0,05=0,3P1+P01(1-P1)
P01=0,0368421
P=(0.9631580,03684210,70,3)

(Sie können die Richtigkeit Ihrer Transtionsmatrix überprüfen, indem Sie sie auf eine hohe Leistung erhöhen. In diesem Fall erledigt 14 die Aufgabe. Jede Zeile des Ergebnisses ergibt die identischen Wahrscheinlichkeiten für den stationären Zustand.)

Beginnen Sie nun in Ihrem Zufallszahlenprogramm mit der zufälligen Auswahl von Zustand 0 oder 1; Dadurch wird ausgewählt, welche Reihe von Sie verwenden. Verwenden Sie dann eine einheitliche Zufallszahl, um den nächsten Zustand zu bestimmen. Diese Nummer ausspucken, ausspülen und nach Bedarf wiederholen.P

Mike Anderson
quelle
Interessante Lösung! Haben Sie vielleicht einen Beispielcode in R? Antone noch?
user333
@Mike Kannst du bitte deinen Account registrieren? Sie sind ein ziemlich aktiver Benutzer, und wir müssen ihn immer wieder manuell zusammenführen. Der Prozess ist ziemlich einfach; Besuchen Sie
Vielen Dank. Wie kann ich anhand der Daten die Markov-Kette (Übergangsmatrix) abschätzen? Gibt es dafür eine R-Funktion?
User333
6

Ich habe einen Riss beim Codieren der Antwort von @Mike Anderson in R gemacht. Ich konnte nicht herausfinden, wie es mit sapply gemacht wird, also habe ich eine Schleife verwendet. Ich habe die Probs leicht verändert, um ein interessanteres Ergebnis zu erzielen, und ich habe 'A' und 'B' verwendet, um die Zustände darzustellen. Lass mich wissen was du denkst.

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: Als Antwort auf Pauls Kommentar hier eine elegantere Formulierung

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

Ich habe den Originalcode geschrieben, als ich gerade R gelernt habe. ;-)

So würden Sie die Übergangsmatrix unter Berücksichtigung der Reihe schätzen:

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

Die Reihenfolge wird gegen meine ursprüngliche Übergangsmatrix ausgetauscht, es werden jedoch die richtigen Wahrscheinlichkeiten ermittelt.

Zach
quelle
Groß! Ich werde so schnell wie möglich ... Sieht gut genug aus ....
user333
Ist es möglich, das Gegenteil zu tun? Angesichts der Reihenschätzung der Matrix?
User333
"Schätzen" ist hier das Schlüsselwort. Die Einträge in der Markov-Kette sind nichts anderes als die bedingten Wahrscheinlichkeiten . Das heißt, Sie können einfach die Reihe scannen und die Proportionen von Nullen, gefolgt von Nullen / Einsen und Einsen, gefolgt von Nullen / Einsen finden, um jede Zeile in der Matrix zu schätzen. Wenn Sie Zeile für Zeile fortfahren, erhalten Sie die üblichen Konfidenzintervalle. Gute Umsetzung. Pr(Xt=ich|Xt-1=j)
Mike Anderson
+1, aber ich habe auch ein paar Kommentare: Eine forSchleife wäre hier etwas sauberer, man kennt die Länge von Series, also einfach benutzen for(i in 2:length(Series)). Dadurch entfällt die Notwendigkeit für i = i + 1. Warum also zuerst probieren Aund dann konvertieren 0,1? Sie könnten direkt Probe 0‚s und 1‘ s.
Paul Hiemstra
2
Im Allgemeinen können Sie es dann in eine neue Funktion 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);
einschließen
1

Hier ist eine Antwort basierend auf dem markovchainPaket, das auf komplexere Abhängigkeitsstrukturen verallgemeinert werden kann.

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)

Dies gibt Ihnen:

Bildbeschreibung hier eingeben

tchakravarty
quelle
0

Ich habe den Überblick über das Papier verloren, in dem dieser Ansatz beschrieben wurde, aber hier ist es.

Zerlegen Sie die Übergangsmatrix in

T=(1-pt)[1001]+pt[p0p0(1-p0)(1-p0)]=(1-pt)ich+ptE

1-ptptp0

ptT11T11=(1-pt)+pt(1-p0)

Eines der nützlichen Merkmale dieser Zerlegung ist, dass sie sich ziemlich einfach auf die Klasse von korrelierten Markov-Modellen in höherdimensionalen Problemen verallgemeinert.

Dave
quelle
Wenn jemand das Papier gesehen hat, das diese Darstellung entwickelt, lassen Sie es mich bitte wissen.
Dave