Test der Signifikanz von Peaks in der Spektraldichte

Wir verwenden manchmal Spektraldichtediagramme, um die Periodizität in Zeitreihen zu analysieren. Normalerweise analysieren wir den Plot durch visuelle Inspektion und versuchen dann, eine Schlussfolgerung über die Periodizität zu ziehen. Aber haben die Statistiker einen Test entwickelt, um zu überprüfen, ob sich Spitzen im Diagramm statistisch von weißem Rauschen unterscheiden? Haben die R-Experten ein Paket für die Spektraldichteanalyse und für diese Art von Tests entwickelt? Großartig, wenn jemand helfen könnte.

Grüße,
P.

Sie sollten sich darüber im Klaren sein, dass das Schätzen von Leistungsspektren mit einem Periodogramm nicht empfohlen wird und in der Tat seit ~ 1896 eine schlechte Praxis ist. Es ist ein inkonsistenter Schätzer für weniger als Millionen von Datenproben (und selbst dann ...) und im Allgemeinen voreingenommen. Das gleiche gilt für die Verwendung von Standardschätzungen für Autokorrelationen (z. B. Bartlett), da es sich um Fourier-Transformationspaare handelt. Vorausgesetzt, Sie verwenden einen konsistenten Schätzer, stehen Ihnen einige Optionen zur Verfügung.

Das Beste davon ist eine Schätzung der Leistungsspektren mit mehreren Fenstern (oder Verjüngung). In diesem Fall können Sie unter Verwendung der Koeffizienten jedes Fensters bei einer interessierenden Frequenz eine harmonische F-Statistik anhand einer Nullhypothese für weißes Rauschen berechnen . Dies ist ein hervorragendes Werkzeug zur Erkennung von Leitungskomponenten im Rauschen und wird dringend empfohlen. Dies ist die Standardeinstellung in der signalverarbeitenden Community für die Erkennung von Periodizitäten im Rauschen unter der Annahme der Stationarität.

Sie können sowohl auf die Multitaper-Methode zur Spektrumsschätzung als auch auf den zugehörigen F-Test über das multitaperPaket in R zugreifen (verfügbar über CRAN). Die mit dem Paket gelieferte Dokumentation sollte ausreichen, um Sie zum Laufen zu bringen. Der F-Test ist eine einfache Option im Funktionsaufruf für spec.mtm.

Die ursprüngliche Referenz, die diese beiden Techniken definiert und die Algorithmen für sie angibt , ist Spectrum Estimation and Harmonic Analysis , DJ Thomson, Proceedings of the IEEE, vol. 70, pg. 1055-1096, 1982.

Hier ist ein Beispiel mit dem mitgelieferten Datensatz multitaper.

require(multitaper);
data(willamette);
resSpec <- spec.mtm(willamette, k=10, nw=5.0, nFFT = "default",
                    centreWithSlepians = TRUE, Ftest = TRUE,
                    jackknife = FALSE, maxAdaptiveIterations = 100,
                    plot = TRUE, na.action = na.fail) 

Die Parameter, die Sie kennen sollten, sind k und nw : Dies sind die Anzahl der Fenster (oben auf 10 eingestellt) und das Zeit-Bandbreiten-Produkt (oben auf 5,0). Sie können diese Werte für die meisten Anwendungen problemlos auf diesen Quasi-Standardwerten belassen. Der Befehl centreWithSlepians entfernt eine robuste Schätzung des Mittelwerts der Zeitreihe mithilfe einer Projektion auf Slepian-Fenster - dies wird auch empfohlen, da das Belassen des Mittelwerts bei niedrigen Frequenzen viel Leistung erzeugt.

Ich würde auch empfehlen, die Spektrum-Ausgabe von 'spec.mtm' auf einer Log-Skala zu zeichnen, da dies die Dinge erheblich aufräumt. Wenn Sie weitere Informationen benötigen, senden Sie einfach eine E-Mail und ich gebe diese gerne weiter.

Wir haben kürzlich in diesem Artikel versucht, dieses Problem durch eine Wavelet-Transformation eines spektralbasierten Tests zu lösen . Im Wesentlichen müssen Sie die Verteilung der Periodogramm-Ordinaten berücksichtigen, ähnlich dem Artikel von Fisher, der in den früheren Antworten erwähnt wurde. Ein anderes Papier von Koen ist dieses . Wir haben kürzlich ein R-Paket hwwntest veröffentlicht .

$f(\omega_k)$

Weitere Einzelheiten zum Test finden Sie in MB Priestley, Spektralanalyse und Zeitreihen , Academic Press, London, 1981, Seite 406.

In R enthält das Paket GeneCycle die Funktion fisher.g.test():

library(GeneCycle)
?fisher.g.test

Hoffe das hilft.