Teste da significância dos picos na densidade espectral

20

Às vezes usamos gráficos de densidade espectral para analisar a periodicidade em séries temporais. Normalmente, analisamos o gráfico por inspeção visual e, em seguida, tentamos tirar uma conclusão sobre a periodicidade. Mas os estatísticos desenvolveram algum teste para verificar se algum pico no gráfico é estatisticamente diferente do ruído branco? Os especialistas da R desenvolveram algum pacote para análise de densidade espectral e para fazer esse tipo de teste? Ótimo se alguém pudesse ajudar.

Atenciosamente,
P.

r time-series hypothesis-testing Pantera
fonte
1
Pressionado por @Wesley, eu apaguei meus pensamentos rápidos sobre funções de autocorrelação e periodograma (pode ser que ele realmente seja um guru de análise de domínio de frequência, mas eu pessoalmente não acho que Bartlett esteja trabalhando com autocorrelações no domínio do tempo), mas ainda acho que meu segunda sugestão sobre bootspecdenspode ser útil.
Dmitrij Celov
Estou baseando minha suposição sobre a resposta das pessoas a 'o que é uma autocorrelação?' nas aparências da literatura, nas quais quase todos os casos em que uma autocorrelação é usada são do padrão, computado no domínio do tempo, autocorrelação de Barlett. E, infelizmente, isso é ruim! :) Agradeço a sugestão de bootspecdensDmitrij; ansioso para verificá-lo.
Wesley Burr

Respostas:

9

Você deve estar ciente de que estimar espectros de potência usando um periodograma não é recomendado e, de fato, é uma prática ruim desde ~ 1896. É um estimador inconsistente para algo menos que milhões de amostras de dados (e mesmo assim ...) e geralmente tendencioso. O mesmo se aplica ao uso de estimativas padrão de autocorrelações (por exemplo, Bartlett), pois são pares de transformadas de Fourier. Desde que você esteja usando um estimador consistente, existem algumas opções disponíveis.

O melhor deles é uma estimativa de várias janelas (ou conicidade) dos espectros de potência. Nesse caso, usando os coeficientes de cada janela com uma frequência de interesse, você pode calcular uma estatística F harmônica contra uma hipótese nula de ruído branco. Essa é uma excelente ferramenta para detecção de componentes de linha em ruídos e é altamente recomendada. É a opção padrão na comunidade de processamento de sinal para detecção de periodicidades no ruído, sob suposição de estacionariedade.

Você pode acessar o método multitaper de estimativa de espectro e o teste F associado através do multitaperpacote em R (disponível via CRAN). A documentação que acompanha o pacote deve ser suficiente para você prosseguir; o teste F é uma opção simples na chamada de função spec.mtm.

A referência original que define essas duas técnicas e fornece os algoritmos para elas é a Estimativa de Espectro e Análise Harmônica , DJ Thomson, Proceedings of the IEEE, vol. 70, pág. 1055-1096, 1982.

Aqui está um exemplo usando o conjunto de dados incluído no multitaperpacote.

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)

Os parâmetros dos quais você deve estar ciente são k e nw : são o número de janelas (definido como 10 acima) e o produto de largura de banda temporal (5.0 acima). Você pode facilmente deixá-los nesses valores quase padrão para a maioria dos aplicativos. O comando centreWithSlepians remove uma estimativa robusta da média da série temporal usando uma projeção nas janelas do Slepian - isso também é recomendado, pois deixar a média produz muita energia nas baixas frequências.

Eu também recomendaria plotar a saída do espectro de 'spec.mtm' em uma escala de log, pois isso limpa significativamente as coisas. Se precisar de mais informações, basta postar e será um prazer fornecer.

Wesley Burr
fonte
A Burr, Silva e Celov - muito obrigado por suas respostas e sugestões interessantes. Estou ansioso para testar esses estimadores. Atenciosamente
Pantera
(+1) nesta noite, pensei cuidadosamente em suas sugestões e decidi que o domínio do tempo é realmente a última coisa (devido ao truncamento do atraso e às propriedades fracas em pequenas amostras) para tentar procurar o comportamento do ciclismo. Pessoalmente, estou preocupado com as suposições para as estatísticas F e as pequenas propriedades de tamanho de amostra do esquema sugerido. Bem, e provavelmente é bom iniciar uma pergunta separada sobre a seleção ideal de janelas, porque existem de fato muitas.
Dmitrij Celov
De fato, existem muitas opções de janela, embora as duas mais comuns sejam as Seqüências Esferoidais Proladas Discretas (ou Slepians ) e as afunilamentos . Se você está procurando a concentração máxima de energia em uma largura de banda local, os Slepianos têm se mostrado ótimos e, de fato, são a saída da forma da equação integral da densidade espectral (consulte o artigo que mencionei para obter detalhes completos). No que diz respeito às estatísticas F, certamente existem alguns problemas com graus de liberdade, mas no geral eles funcionam muito bem, com ~ 2k-2 dof disponível.
Wesley Burr
Periodograma suavizado também usa cone, permite FFT, o livro de David Stoffer ensina como calcular intervalos de confiança também. Este multitaperpacote parece ter empregado técnicas mais avançadas para diminuir e calcular o intervalo de confiança. Mas acho que a ideia era a mesma, de acordo com David Stoffer. Esta é a única coisa que eu poderia pensar que ensinar peridogoram de baunilha ainda faz sentido hoje.
stucash
ok, então você é um dos autores deste pacote e usou algumas palavras muito fortes contra o periodograma. Espero que um dia você possa voltar com mais evidências. Prós e contras comuns do Periodograma são bem conhecidos, como sua variação explosiva, motivo pelo qual não é um bom estimador consistente do espectro, mas o periodograma suavizado não é tão ruim assim, como você afirmou aqui.
stucash
3

Tentamos resolver esse problema com uma transformação wavelet de um teste espectral recentemente neste artigo . Essencialmente, você precisa considerar a distribuição de ordenadas do periodograma, da mesma forma que o artigo de Fisher, mencionado nas respostas anteriores. Outro artigo de Koen é este . Publicamos recentemente um pacote R hwwntest .

Delyan Savchev
fonte
Savchev, muito obrigado pelo seu comentário e referências. Estou ansioso para testar seu pacote R.
Pantera20 de
2

f(ωk)

Você pode obter mais detalhes sobre o teste em MB Priestley, Análise espectral e séries temporais , Academic Press, Londres, 1981, página 406.

Em R, o pacote GeneCycle contém a função fisher.g.test():

library(GeneCycle)
?fisher.g.test

Espero que isto ajude.

Washington S. Silva
fonte
isso é ótimo, mas o teste g do pacote se baseia em sua própria função de periodograma, que tem opções muito limitadas para calcular espectros de potência ...
stucash