Deteksi periode dari deret waktu umum

53

Posting ini adalah kelanjutan dari posting lain yang terkait dengan metode umum untuk deteksi outlier dalam deret waktu . Pada dasarnya, pada titik ini saya tertarik pada cara yang kuat untuk menemukan periodisitas / musim dari rangkaian waktu generik yang dipengaruhi oleh banyak kebisingan. Dari sudut pandang pengembang, saya ingin antarmuka yang sederhana seperti:

unsigned int discover_period(vector<double> v);

Di mana varray berisi sampel, dan nilai kembali adalah periode sinyal. Poin utamanya adalah, sekali lagi, saya tidak bisa membuat asumsi mengenai sinyal yang dianalisis. Saya sudah mencoba pendekatan yang didasarkan pada autokorelasi sinyal (mendeteksi puncak suatu korelasiogram), tetapi tidak sekuat yang saya inginkan.

gianluca
sumber
1
Sudahkah Anda mencoba xts :: periodisitas?
Fabrício

Jawaban:

49

Jika Anda benar-benar tidak tahu apa periodisitasnya, mungkin pendekatan terbaik adalah menemukan frekuensi yang sesuai dengan maksimum kepadatan spektral. Namun, spektrum pada frekuensi rendah akan dipengaruhi oleh tren, jadi Anda harus membatalkan seri terlebih dahulu. Fungsi R berikut harus melakukan pekerjaan untuk sebagian besar seri. Ini jauh dari sempurna, tetapi saya telah mengujinya pada beberapa lusin contoh dan tampaknya berfungsi ok. Ini akan mengembalikan 1 untuk data yang tidak memiliki periodisitas kuat, dan panjang periode sebaliknya.

Pembaruan: Versi 2 fungsi. Ini jauh lebih cepat dan tampaknya lebih kuat.

find.freq <- function(x)
{
    n <- length(x)
    spec <- spec.ar(c(x),plot=FALSE)
    if(max(spec$spec)>10) # Arbitrary threshold chosen by trial and error.
    {
        period <- round(1/spec$freq[which.max(spec$spec)])
        if(period==Inf) # Find next local maximum
        {
            j <- which(diff(spec$spec)>0)
            if(length(j)>0)
            {
                nextmax <- j[1] + which.max(spec$spec[j[1]:500])
                period <- round(1/spec$freq[nextmax])
            }
            else
                period <- 1
        }
    }
    else
        period <- 1
    return(period)
}
Rob Hyndman
sumber
Terima kasih. Sekali lagi, saya akan mencoba pendekatan ini sesegera mungkin dan akan menulis di sini hasil akhirnya.
gianluca
2
Gagasan Anda cukup bagus, tetapi dalam kasus saya, gagal mendeteksi periodisitas rangkaian waktu yang sangat sederhana (dan tidak terlalu berisik) seperti dl.dropbox.com/u/540394/chart.png . Dengan pendekatan "empiris" (berdasarkan autokorelasi), algoritma sederhana yang saya tulis mengembalikan periode tepat 1008 (memiliki sampel setiap 10 menit, ini berarti 1008/24/6 = 7, jadi periodisitas mingguan). Masalah utama saya adalah: 1) Terlalu lambat untuk berkumpul (ini membutuhkan banyak data historis) dan saya memerlukan pendekatan online reaktif; 2) Sangat tidak efisien dari sudut pandang penggunaan memori; 3) Tidak kuat sama sekali;
gianluca
Terima kasih. Sayangnya, ini masih tidak berfungsi seperti yang saya harapkan. Untuk rangkaian waktu yang sama dari komentar sebelumnya, ia mengembalikan 166, yang hanya sebagian benar (dari sudut pandang saya, periode mingguan yang jelas lebih menarik). Dan menggunakan deret waktu yang sangat bising, seperti ini dl.dropbox.com/u/540394/chart2.png (analisis jendela penerima TCP), fungsinya mengembalikan 10, sementara saya mengharapkan 1 (saya tidak dapat melihat dengan jelas periodisitas). BTW Saya tahu bahwa akan sangat sulit untuk menemukan apa yang saya cari, karena saya berurusan dengan sinyal yang terlalu berbeda.
gianluca
166 bukan perkiraan buruk 168. Jika Anda tahu data diamati setiap jam dengan pola mingguan, lalu mengapa memperkirakan frekuensi sama sekali?
Rob Hyndman
5
Versi yang ditingkatkan dalam paket perkiraan sebagaifindfrequency
Rob Hyndman
10

Jika Anda mengharapkan proses menjadi diam - periodisitas / musim tidak akan berubah seiring waktu - maka sesuatu seperti periodogram Chi-square (lihat misalnya Sokolove dan Bushell, 1978) mungkin merupakan pilihan yang baik. Ini umumnya digunakan dalam analisis data sirkadian yang dapat memiliki jumlah suara yang sangat besar di dalamnya, tetapi diharapkan memiliki periodikitas yang sangat stabil.

Pendekatan ini tidak membuat asumsi tentang bentuk gelombang (selain dari itu konsisten dari siklus ke siklus), tetapi memang mensyaratkan bahwa setiap kebisingan menjadi rata-rata konstan dan tidak berkorelasi dengan sinyal.

chisq.pd <- function(x, min.period, max.period, alpha) {
N <- length(x)
variances = NULL
periods = seq(min.period, max.period)
rowlist = NULL
for(lc in periods){
    ncol = lc
    nrow = floor(N/ncol)
    rowlist = c(rowlist, nrow)
    x.trunc = x[1:(ncol*nrow)]
    x.reshape = t(array(x.trunc, c(ncol, nrow)))
    variances = c(variances, var(colMeans(x.reshape)))
}
Qp = (rowlist * periods * variances) / var(x)
df = periods - 1
pvals = 1-pchisq(Qp, df)
pass.periods = periods[pvals<alpha]
pass.pvals = pvals[pvals<alpha]
#return(cbind(pass.periods, pass.pvals))
return(cbind(periods[pvals==min(pvals)], pvals[pvals==min(pvals)]))
}

x = cos( (2*pi/37) * (1:1000))+rnorm(1000)
chisq.pd(x, 2, 72, .05)

Dua baris terakhir hanyalah sebuah contoh, yang menunjukkan bahwa ia dapat mengidentifikasi periode fungsi trigonometri murni, bahkan dengan banyak noise tambahan.

Seperti yang ditulis, argumen terakhir ( alpha) dalam panggilan itu berlebihan, fungsinya hanya mengembalikan periode 'terbaik' yang dapat ditemukan; batalkan komentar pertama returndan komentar kedua untuk mengembalikan daftar semua periode yang signifikan di level tersebut alpha.

Fungsi ini tidak melakukan pengecekan kewarasan apa pun untuk memastikan bahwa Anda telah memasukkan periode yang dapat diidentifikasi, juga tidak (dapat) bekerja dengan periode fraksional, juga tidak ada semacam kontrol perbandingan ganda yang dibangun jika Anda memutuskan untuk lihat beberapa periode. Tetapi selain itu harus cukup kuat.

Kaya
sumber
Terlihat menarik tapi saya tidak mengerti hasilnya, tidak memberi tahu saya di mana periode dimulai, dan sebagian besar nilai dari 1.
Herman Toothrot
3

Anda mungkin ingin mendefinisikan apa yang Anda inginkan dengan lebih jelas (untuk diri sendiri, jika tidak di sini). Jika yang Anda cari adalah periode diam paling signifikan secara statistik yang terkandung dalam data bising Anda, pada dasarnya ada dua rute yang harus diambil:

1) menghitung estimasi autokorelasi yang kuat, dan mengambil koefisien maksimum
2) menghitung estimasi kepadatan spektral daya yang kuat, dan mengambil maksimum spektrum

Masalah dengan # 2 adalah untuk deret waktu yang berisik, Anda akan mendapatkan sejumlah besar daya dalam frekuensi rendah, sehingga sulit untuk dibedakan. Ada beberapa teknik untuk menyelesaikan masalah ini (yaitu pra-pemutihan, lalu perkirakan PSD), tetapi jika periode sebenarnya dari data Anda cukup lama, deteksi otomatis akan rapuh.

Taruhan terbaik Anda mungkin untuk menerapkan rutin autokorelasi yang kuat seperti yang dapat ditemukan di bab 8.6, 8.7 di Robust Statistics - Theory and Methods oleh Maronna, Martin dan Yohai. Mencari Google untuk "durbin-levinson yang kuat" juga akan menghasilkan beberapa hasil.

Jika Anda hanya mencari jawaban sederhana, saya tidak yakin jawabannya ada. Deteksi periode dalam deret waktu bisa rumit, dan meminta rutin otomatis yang dapat melakukan sihir mungkin terlalu banyak.

Wesley Burr
sumber
Terima kasih atas informasi berharga Anda, saya akan melihat buku itu dengan pasti.
gianluca
3

Anda bisa menggunakan Hilbert Transformation from the DSP theory untuk mengukur frekuensi sesaat dari data Anda. Situs http://ta-lib.org/ memiliki kode sumber terbuka untuk mengukur periode siklus dominan dari data keuangan; fungsi yang relevan disebut HT_DCPERIOD; Anda mungkin dapat menggunakan ini atau menyesuaikan kode dengan tujuan Anda.

babelproofreader
sumber
3

Pendekatan yang berbeda bisa berupa Dekomposisi Mode Empiris. Paket R disebut EMD yang dikembangkan oleh penemu metode:

require(EMD)
ndata <- 3000  
tt2 <- seq(0, 9, length = ndata)  
xt2 <- sin(pi * tt2) + sin(2* pi * tt2) + sin(6 * pi * tt2) + 0.5 * tt2  
try <- emd(xt2, tt2, boundary = "wave")  
### Ploting the IMF's  
par(mfrow = c(try$nimf + 1, 1), mar=c(2,1,2,1))  
rangeimf <- range(try$imf)  
for(i in 1:try$nimf) {  
plot(tt2, try$imf[,i], type="l", xlab="", ylab="", ylim=rangeimf, main=paste(i, "-th IMF", sep="")); abline(h=0)  
}  
plot(tt2, try$residue, xlab="", ylab="", main="residue", type="l", axes=FALSE); box()

Metode ini diberi merek 'Empiris' untuk alasan yang baik dan ada risiko bahwa Fungsi Mode Intrinsik (komponen aditif individu) ikut campur. Di sisi lain metode ini sangat intuitif dan dapat membantu untuk inspeksi visual yang cepat dari siklus.

Fabrizio Maccallini
sumber
0

Mengacu pada pos Rob Hyndman di atas https://stats.stackexchange.com/a/1214/70282

Fungsi find.freq bekerja dengan sangat baik. Pada set data harian yang saya gunakan, frekuensi frekuensinya menjadi 7.

Ketika saya mencobanya hanya pada hari-hari minggu, frekuensinya adalah 23, yang sangat dekat dengan 21,42857 = 29,6 * 5/7 yang merupakan jumlah rata-rata hari kerja dalam sebulan. (Atau sebaliknya 23 * 7/5 adalah 32.)

Melihat kembali data harian saya, saya bereksperimen dengan firasat mengambil periode pertama, rata-rata dengan itu dan kemudian menemukan periode berikutnya, dll. Lihat di bawah:

find.freq.all = fungsi (x) {  
  f = find.freq (x);
  freqs = c (f);  
  while (f> 1) {
    mulai = 1; #juga coba mulai = f;
    x = period.apply (x, seq (mulai, panjang (x), f), rata-rata); 
    f = find.freq (x);
    freqs = c (freqs, f);
  }
  if (length (freqs) == 1) {return (freqs); }
  untuk (i in 2: length (freqs)) {
    freqs [i] = freqs [i] * freqs [i-1];
  }
  freqs [1: (length (freqs) -1)];
}
find.freq.all (dailyts) #menggunakan data harian

Di atas memberi (7,28) atau (7,35) tergantung pada apakah seq dimulai dengan 1 atau f. (Lihat komentar di atas.)

Yang akan menyiratkan bahwa periode musiman untuk msts (...) harus (7,28) atau (7,35).

Logika muncul sensitif terhadap kondisi awal mengingat sensitivitas parameter algoritma. Rata-rata 28 dan 35 adalah 31,5 yang dekat dengan panjang rata-rata sebulan.

Saya kira saya menemukan kembali roda, apa nama algoritma ini? Apakah ada implementasi yang lebih baik di R di suatu tempat?

Kemudian, saya menjalankan kode di atas dalam mencoba semua mulai dari 1 sampai 7 dan saya mendapat 35,35,28,28,28,28,28 untuk periode kedua. Rata-rata berhasil hingga 30 yang merupakan jumlah rata-rata hari dalam sebulan. Menarik...

Ada pemikiran atau komentar?

Chris
sumber
0

Seseorang juga dapat menggunakan tes Ljung-Box untuk mengetahui perbedaan musim mana yang mencapai stasioneritas terbaik. Saya sedang mengerjakan subjek yang berbeda dan saya menggunakan ini sebenarnya untuk tujuan yang sama. Coba periode yang berbeda seperti 3 hingga 24 untuk data bulanan. Dan uji masing-masing dengan Ljung-Box dan simpan hasil Chi-Square. Dan pilih periode dengan nilai chi-square terendah.

Berikut adalah kode sederhana untuk melakukannya.

minval0 <- 5000 #assign a big number to be sure Chi values are smaller
minindex0 <- 0
periyot <- 0

for (i in 3:24) { #find optimum period by Qtests over original data

        d0D1 <- diff(a, lag=i)

        #store results
        Qtest_d0D1[[i]] <- Box.test(d0D1, lag=20, type = "Ljung-Box")

        #store Chi-Square statistics
        sira0[i] <- Qtest_d0D1[[i]][1]
}
#turn list to a data frame, then matrix
datam0 <- data.frame(matrix(unlist(sira0), nrow=length(Qtest_d0D1)-2, byrow = T))
datamtrx0 <- as.matrix(datam0[])
#get min value's index
minindex0 <- which(datamtrx0 == min(datamtrx0), arr.ind = F)
periyot <- minindex0 + 2
ali
sumber