Bagaimana cara menghasilkan data deret waktu biner acak yang dihubungkan secara acak?

15

Bagaimana saya bisa menghasilkan deret waktu biner sehingga:

  1. Probabilitas rata-rata mengamati 1 ditentukan (misalnya 5%);
  2. Probabilitas kondisional untuk mengamati 1 pada waktu diberi nilai pada t - 1 (katakanlah 30% jika nilai t - 1 adalah 1)?tt1t1
pengguna333
sumber

Jawaban:

17

Gunakan rantai Markov dua negara.

Jika keadaan disebut 0 dan 1, maka rantai dapat diwakili oleh matriks 2x2 memberikan probabilitas transisi antara keadaan, di mana P i j adalah probabilitas untuk berpindah dari keadaan i ke keadaan j . Dalam matriks ini, setiap baris harus berjumlah 1,0.PPijij

Dari pernyataan 2, kita memiliki , dan konservasi sederhana kemudian mengatakan P 10 = 0,7 .P11=0.3P10=0.7

Dari pernyataan 1, Anda ingin probabilitas jangka panjang (juga disebut keseimbangan atau kondisi-mapan) menjadi . Ini mengatakan P 1 = 0,05 = 0,3 P 1 + P 01 ( 1 - P 1 ) Memecahkan memberikan P 01 = 0,0368421 dan matriks transisi 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)

(Anda dapat memeriksa matriks transtion Anda dengan benar dengan menaikkannya ke daya tinggi - dalam hal ini 14 melakukan pekerjaan - setiap baris hasilnya memberikan probabilitas kondisi steady state yang identik)

P

Mike Anderson
sumber
Solusi menarik! Apakah Anda mungkin memiliki beberapa kode sampel dalam R? Antone lain?
user333
@ Mike Bisakah Anda mendaftarkan akun Anda? Anda adalah pengguna yang cukup aktif dan kami harus menggabungkannya secara manual berulang kali. Prosesnya cukup mudah; cukup kunjungi stats.stackexchange.com/login
Terima kasih. Bagaimana saya bisa memperkirakan rantai Markov (matriks transisi) yang diberikan data? Apakah ada fungsi R untuk melakukan itu?
user333
6

Saya mengambil celah pada coding @ Mike Anderson jawaban dalam R. Saya tidak tahu bagaimana melakukannya menggunakan sapply, jadi saya menggunakan loop. Saya sedikit mengubah probs untuk mendapatkan hasil yang lebih menarik, dan saya menggunakan 'A' dan 'B' untuk mewakili negara bagian. Biarkan aku tahu apa yang kamu pikirkan.

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: Menanggapi komentar Paul, berikut ini formulasi yang lebih elegan

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

Saya menulis kode asli ketika saya baru belajar R, jadi potong sedikit kendur. ;-)

Inilah cara Anda memperkirakan matriks transisi, mengingat seri:

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

Urutan ditukar vs matriks transisi asli saya, tetapi mendapat probabilitas yang tepat.

Zach
sumber
Bagus! Saya akan sesegera mungkin ... Terlihat cukup baik ....
user333
Apakah mungkin untuk melakukan invers? Diberikan seri estimasi matriks?
user333
Pr(Xt=i|Xt1=j)
+1, tapi saya juga punya beberapa komentar: Sebuah forloop akan menjadi sedikit lebih bersih di sini, Anda tahu panjangnya Series, jadi gunakan saja for(i in 2:length(Series)). Ini menghilangkan kebutuhan untuk i = i + 1. Juga, mengapa sampel pertama A, lalu dikonversi menjadi 0,1? Anda dapat langsung mencicipi 0dan 1.
Paul Hiemstra
2
Secara umum, Anda dapat membungkusnya dengan fungsi baru 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);untuk memungkinkan autokorelasi lag 1 yang ditentukan sebelumnya
Tom Wenseleers
1

Berikut adalah jawaban berdasarkan markovchainpaket yang dapat digeneralisasikan ke struktur ketergantungan yang lebih kompleks.

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)

Ini memberi Anda:

masukkan deskripsi gambar di sini

tchakravarty
sumber
0

Saya kehilangan jejak makalah di mana pendekatan ini dijelaskan, tapi begini saja.

Menguraikan matriks transisi menjadi

T=(1-halt)[1001]+halt[hal0hal0(1-hal0)(1-hal0)]=(1-halt)saya+haltE

yang, secara intuitif, sesuai dengan gagasan bahwa ada beberapa kemungkinan 1-halt bahwa sistem tetap dalam keadaan yang sama, dan probabilitas halt bahwa keadaan menjadi acak, di mana acak berarti membuat undian independen dari distribusi kesetimbangan untuk keadaan berikutnya (hal0 adalah probabilitas keseimbangan untuk berada dalam keadaan pertama).

Perhatikan bahwa dari data yang Anda tentukan perlu Anda pecahkan halt dari yang ditentukan T11 melalui T11=(1-halt)+halt(1-hal0).

Salah satu fitur yang berguna dari dekomposisi ini adalah bahwa itu cukup mudah digeneralisasikan ke kelas model Markov berkorelasi dalam masalah dimensi yang lebih tinggi.

Dave
sumber
Jika ada yang melihat makalah yang mengembangkan representasi ini, beri tahu saya.
Dave