Wykrywanie okresu ogólnych szeregów czasowych

Ten post jest kontynuacją kolejnego postu związanego z ogólną metodą wykrywania wartości odstających w szeregach czasowych . Zasadniczo w tym momencie interesuje mnie solidny sposób odkrywania okresowości / sezonowości ogólnych szeregów czasowych dotkniętych dużym hałasem. Z punktu widzenia programisty chciałbym prosty interfejs, taki jak:

unsigned int discover_period(vector<double> v);

Gdzie vjest tablica zawierająca próbki, a zwracana wartość to okres sygnału. Chodzi o to, że znowu nie mogę przyjąć żadnego założenia dotyczącego analizowanego sygnału. Próbowałem już podejścia opartego na autokorelacji sygnału (wykrywanie szczytów korelogramu), ale nie jest ono solidne, jak bym chciał.

Jeśli naprawdę nie masz pojęcia, co to jest okresowość, prawdopodobnie najlepszym rozwiązaniem jest znalezienie częstotliwości odpowiadającej maksymalnej gęstości widmowej. Jednak na spektrum przy niskich częstotliwościach będzie miał wpływ trend, więc najpierw musisz odrzucić serię. Następująca funkcja R powinna wykonać zadanie dla większości serii. Jest daleki od ideału, ale przetestowałem go na kilkudziesięciu przykładach i wydaje się, że działa dobrze. Zwróci 1 dla danych, które nie mają silnej okresowości, a długość okresu w przeciwnym razie.

Aktualizacja: Wersja 2 funkcji. Jest to znacznie szybsze i wydaje się bardziej niezawodne.

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

Jeśli oczekujesz, że proces będzie stacjonarny - okresowość / sezonowość nie zmieni się w czasie - wtedy coś w rodzaju okresogramu chi-kwadrat (patrz np. Sokolove i Bushell, 1978) może być dobrym wyborem. Jest powszechnie stosowany w analizie danych okołodobowych, które mogą mieć bardzo duże ilości hałasu, ale oczekuje się, że będą miały bardzo stabilne okresy.

Podejście to nie zakłada założenia kształtu fali (poza tym, że jest on spójny między cyklami), ale wymaga, aby każdy szum miał stałą średnią i nie był skorelowany z sygnałem.

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)

Ostatnie dwa wiersze są tylko przykładem, pokazującym, że może on identyfikować okres czystej funkcji trygonometrycznej, nawet przy dużej ilości szumów addytywnych.

Jak napisano, ostatni argument ( alpha) w wywołaniu jest zbędny, funkcja po prostu zwraca „najlepszy” okres, jaki można znaleźć; odkomentuj pierwsze returnstwierdzenie i skomentuj drugie, aby zwrócić listę wszystkich istotnych okresów na poziomie alpha.

Ta funkcja nie dokonuje żadnego sprawdzenia poprawności poczytalności, aby upewnić się, że wprowadziłeś identyfikowalne okresy, ani nie (może to) działać z okresami ułamkowymi, ani nie ma wbudowanej kontroli wielokrotnego porównania, jeśli zdecydujesz się spójrz na wiele okresów. Ale poza tym powinien być dość solidny.

Możesz lepiej zdefiniować, co chcesz (dla siebie, jeśli nie tutaj). Jeśli to, czego szukasz, jest najistotniejszym statystycznie okresem stacjonarnym zawartym w zaszumionych danych, istnieją zasadniczo dwie trasy:

1) obliczyć solidne oszacowanie autokorelacji i przyjąć maksymalny współczynnik
2) obliczyć solidne oszacowanie gęstości widmowej mocy i przyjąć maksimum widma

Problem z numerem 2 polega na tym, że dla każdej hałaśliwej serii czasowej otrzymasz dużą moc na niskich częstotliwościach, co utrudnia rozróżnienie. Istnieją pewne techniki rozwiązania tego problemu (np. Wstępne wybielanie, a następnie oszacowanie PSD), ale jeśli prawdziwy okres na podstawie danych jest wystarczająco długi, automatyczne wykrywanie będzie trudne.

Najlepszym rozwiązaniem jest prawdopodobnie wdrożenie solidnej procedury autokorelacji, takiej jak opisana w rozdziale 8.6, 8.7 w Solidnej statystyki - teorii i metod autorstwa Maronny, Martina i Yohai. Wyszukanie w Google hasła „solidny durbin-levinson” również przyniesie pewne rezultaty.

Jeśli szukasz prostej odpowiedzi, nie jestem pewien, czy istnieje. Wykrywanie okresu w szeregach czasowych może być skomplikowane, a prośba o zautomatyzowaną procedurę, która może wykonywać magię, może być zbyt duża.

Możesz użyć transformacji Hilberta z teorii DSP do pomiaru chwilowej częstotliwości twoich danych. Witryna http://ta-lib.org/ ma otwarty kod źródłowy do pomiaru dominującego okresu cyklu danych finansowych; odpowiednia funkcja nosi nazwę HT_DCPERIOD; możesz to wykorzystać lub dostosować kod do swoich celów.

Innym podejściem może być rozkład w trybie empirycznym. Pakiet R nazywa się EMD opracowany przez wynalazcę metody:

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

Metodę tę nazwano „Empiryczną” nie bez powodu i istnieje ryzyko, że funkcje trybu wewnętrznego (poszczególne składniki dodatków) zostaną pomieszane. Z drugiej strony metoda jest bardzo intuicyjna i może być pomocna w szybkiej wizualnej kontroli cykliczności.

W nawiązaniu do postu Roba Hyndmana powyżej https://stats.stackexchange.com/a/1214/70282

Funkcja find.freq działa doskonale. W codziennym zestawie danych, którego używam, poprawnie wyliczyła częstotliwość na 7.

Kiedy wypróbowałem to tylko w dni robocze, wspomniałem, że częstotliwość wynosi 23, co jest niezwykle zbliżone do 21,42857 = 29,6 * 5/7, co jest średnią liczbą dni roboczych w miesiącu. (Lub odwrotnie: 23 * 7/5 to 32.)

Patrząc wstecz na moje codzienne dane, eksperymentowałem z przeczuciem, biorąc pierwszy okres, uśredniając go, a następnie znajdując następny okres itp. Zobacz poniżej:

find.freq.all = funkcja (x) {  
  f = find.freq (x);
  freqs = c (f);  
  podczas gdy (f> 1) {
    start = 1; # także spróbuj start = f;
    x = period.apply (x, seq (start, length (x), f), mean); 
    f = find.freq (x);
    freqs = c (freqs, f);
  }
  if (length (freqs) == 1) {return (freqs); }
  dla (i in 2: length (freqs)) {
    freqs [i] = freqs [i] * freqs [i-1];
  }
  freqs [1: (length (freqs) -1)];
}
find.freq.all (dailyts) # wykorzystując codzienne dane

Powyższe daje (7,28) lub (7,35) w zależności od tego, czy sekwencja zaczyna się od 1 czy f. (Patrz komentarz powyżej.)

Co oznaczałoby, że okresy sezonowe dla msts (...) powinny wynosić (7,28) lub (7,35).

Logika wydaje się wrażliwa na warunki początkowe, biorąc pod uwagę czułość parametrów algorytmu. Średnia z 28 i 35 wynosi 31,5, co jest zbliżone do średniej długości miesiąca.

Podejrzewam, że wymyśliłem koło na nowo, jak nazywa się ten algorytm? Czy jest gdzieś lepsza implementacja w R?

Później uruchomiłem powyższy kod, próbując wszystkich początków od 1 do 7, i otrzymałem 35.325.28,28,28,28 za drugi okres. Średnia oblicza się do 30, co jest średnią liczbą dni w miesiącu. Ciekawy...

Wszelkie myśli lub komentarze?

Można również użyć testu Ljunga-Boxa, aby dowiedzieć się, która różnica sezonowa osiąga najlepszą stacjonarność. Pracowałem nad innym tematem i wykorzystałem to w rzeczywistości do tych samych celów. Wypróbuj różne okresy, np. Od 3 do 24, aby uzyskać dane miesięczne. Przetestuj każdy z nich przez Ljung-Box i zapisz wyniki Chi-Square. I wybierz okres o najniższej wartości chi-kwadrat.

Oto prosty kod, aby to zrobić.

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