Detecção de período de uma série temporal genérica

53

Este post é a continuação de outro post relacionado a um método genérico para detecção de outlier em séries temporais . Basicamente, neste ponto, estou interessado em uma maneira robusta de descobrir a periodicidade / sazonalidade de uma série temporal genérica afetada por muito ruído. Do ponto de vista do desenvolvedor, eu gostaria de uma interface simples como:

unsigned int discover_period(vector<double> v);

Onde vestá o array que contém as amostras e o valor de retorno é o período do sinal. O ponto principal é que, novamente, não posso fazer nenhuma suposição sobre o sinal analisado. Eu já tentei uma abordagem baseada na autocorrelação de sinal (detectando os picos de um correlograma), mas não é robusta como eu gostaria.

gianluca
fonte
11
Você já experimentou xts :: periodicity?
Fabrício

Respostas:

49

Se você realmente não tem idéia de qual é a periodicidade, provavelmente a melhor abordagem é encontrar a frequência correspondente ao máximo da densidade espectral. No entanto, o espectro em baixas frequências será afetado pela tendência; portanto, você deve prejudicar a série primeiro. A seguinte função R deve fazer o trabalho para a maioria das séries. Está longe de ser perfeito, mas já o testei em algumas dezenas de exemplos e parece funcionar bem. Ele retornará 1 para dados que não possuem periodicidade forte e a duração do período caso contrário.

Atualização: versão 2 da função. Isso é muito mais rápido e parece ser mais robusto.

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
fonte
Obrigado. Mais uma vez, tentarei essa abordagem o mais rápido possível e escreverei aqui os resultados finais.
Gianluca
2
Sua ideia é bastante boa, mas, no meu caso, falha em detectar a periodicidade de uma série cronológica realmente simples (e não tão barulhenta) como dl.dropbox.com/u/540394/chart.png . Com minha abordagem "empírica" ​​(baseada na autocorrelação), o algoritmo simples que escrevi retorna um período exato de 1008 (com uma amostra a cada 10 minutos, isso significa 1008/24/6 = 7, portanto, uma periodicidade semanal). Meus principais problemas são: 1) É muito lento para convergir (requer muitos dados históricos) e eu preciso de uma abordagem on-line reativa; 2) É ineficiente como o inferno, do ponto de vista do uso da memória; 3) Não é robusto;
Gianluca
Obrigado. Infelizmente, isso ainda não funciona como eu esperaria. Nas mesmas séries temporais do comentário anterior, retorna 166, o que é parcialmente certo (do meu ponto de vista, o período semanal evidente é mais interessante). E usando uma série temporal muito barulhenta, como esta dl.dropbox.com/u/540394/chart2.png (uma análise da janela do receptor TCP), a função retorna 10, enquanto eu esperaria 1 (não vejo nada óbvio) periodicidade). BTW, eu sei que será realmente difícil encontrar o que estou procurando, pois estou lidando com sinais muito diferentes.
Gianluca
166 não é uma estimativa ruim de 168. Se você sabe que os dados são observados a cada hora com um padrão semanal, por que estimar a frequência?
Rob Hyndman
5
Uma versão melhorada está no pacote de previsões comofindfrequency
Rob Hyndman
10

Se você espera que o processo seja estacionário - a periodicidade / sazonalidade não mudará ao longo do tempo -, algo como um periodograma do qui-quadrado (veja, por exemplo, Sokolove e Bushell, 1978) pode ser uma boa escolha. É comumente usado na análise de dados circadianos que podem ter quantidades extremamente grandes de ruído, mas espera-se que tenha periodicidades muito estáveis.

Essa abordagem não pressupõe a forma da forma de onda (além de ser consistente de ciclo para ciclo), mas exige que qualquer ruído seja de média constante e não correlacionado com o sinal.

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)

As duas últimas linhas são apenas um exemplo, mostrando que ele pode identificar o período de uma função trigonométrica pura, mesmo com muito ruído aditivo.

Como está escrito, o último argumento ( alpha) na chamada é supérfluo; a função simplesmente retorna o 'melhor' período que pode encontrar; remova o comentário da primeira returndeclaração e comente a segunda para que ela retorne uma lista de todos os períodos significativos no nível alpha.

Essa função não realiza nenhum tipo de verificação de sanidade para garantir que você tenha colocado períodos identificáveis, nem (pode) funcionar com períodos fracionários, nem existe algum tipo de controle de comparação múltiplo incorporado, se você decidir olhe para vários períodos. Mas fora isso, deve ser razoavelmente robusto.

Rico
fonte
Parece interessante, mas eu não entendo a saída, não me diga onde inicia o período, ea maioria dos pvalues de 1.
Herman Toothrot
3

Você pode definir o que deseja de forma mais clara (para você mesmo, se não estiver aqui). Se o que você está procurando é o período estacionário estatisticamente mais significativo contido nos seus dados barulhentos, há essencialmente duas rotas a seguir:

1) calcular uma estimativa robusta de autocorrelação e obter o coeficiente máximo
2) calcular uma estimativa robusta da densidade espectral de potência e obter o máximo do espectro

O problema do número 2 é que, para qualquer série temporal barulhenta, você obtém uma grande quantidade de energia em baixas frequências, dificultando a distinção. Existem algumas técnicas para resolver esse problema (ou seja, pré-branquear e estimar o PSD), mas se o período real dos dados for longo o suficiente, a detecção automática será duvidosa.

Sua melhor aposta é provavelmente implementar uma rotina robusta de autocorrelação, como pode ser encontrada no capítulo 8.6, 8.7 em Estatísticas robustas - teoria e métodos de Maronna, Martin e Yohai. Procurar no Google por "durbin-levinson robusto" também produzirá alguns resultados.

Se você está apenas procurando uma resposta simples, não tenho certeza de que exista. A detecção de períodos em séries temporais pode ser complicada, e solicitar uma rotina automatizada que possa executar mágica pode ser demais.

Wesley Burr
fonte
Obrigado por suas informações preciosas, vou ver esse livro com certeza.
Gianluca
3

Você pode usar a transformação de Hilbert da teoria DSP para medir a frequência instantânea de seus dados. O site http://ta-lib.org/ possui código-fonte aberto para medir o período dominante do ciclo de dados financeiros; a função relevante é chamada HT_DCPERIOD; você pode usar isso ou adaptar o código aos seus propósitos.

babelproofreader
fonte
3

Uma abordagem diferente poderia ser a decomposição do modo empírico. O pacote R é chamado EMD desenvolvido pelo inventor do método:

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

O método foi marcado como 'Empírico' por um bom motivo e existe o risco de que as Funções do Modo Intrínseco (os componentes aditivos individuais) se misturem. Por outro lado, o método é muito intuitivo e pode ser útil para uma inspeção visual rápida da ciclicidade.

Fabrizio Maccallini
fonte
0

Em referência à postagem de Rob Hyndman acima https://stats.stackexchange.com/a/1214/70282

A função find.freq funciona de maneira brilhante. No conjunto de dados diário que estou usando, calculou corretamente a frequência como 7.

Quando o experimentei apenas nos dias da semana, mencionou que a frequência é 23, o que é notavelmente próximo de 21,42857 = 29,6 * 5/7, que é o número médio de dias úteis em um mês. (Ou, inversamente, 23 * 7/5 é 32.)

Revendo meus dados diários, experimentei um palpite de tirar o primeiro período, calcular a média desse valor e, em seguida, encontrar o próximo período, etc. Veja abaixo:

find.freq.all = function (x) {  
  f = find.freq (x);
  freqs = c (f);  
  while (f> 1) {
    start = 1; #também tente iniciar = f;
    x = período.aplicação (x, seq (início, comprimento (x), f), média); 
    f = find.freq (x);
    freqs = c (freqs, f);
  }
  if (comprimento (freqs) == 1) {return (freqs); }
  para (i em 2: comprimento (freqs)) {
    freqs [i] = freqs [i] * freqs [i-1];
  }
  freqs [1: (comprimento (freqs) -1)];
}
find.freq.all (dailyts) # usando dados diários

O exemplo acima fornece (7,28) ou (7,35), dependendo de seq começar com 1 ou f. (Veja o comentário acima.)

O que implicaria que os períodos sazonais para msts (...) fossem (7,28) ou (7,35).

A lógica parece sensível às condições iniciais, dada a sensibilidade dos parâmetros do algoritmo. A média de 28 e 35 é 31,5, o que está próximo da duração média de um mês.

Eu suspeito que reinventei a roda, qual é o nome desse algoritmo? Existe uma melhor implementação em R em algum lugar?

Mais tarde, executei o código acima ao tentar todas as partidas de 1 a 7 e obtive 35,35,28,28,28,28,28 pelo segundo período. A média chega a 30, que é o número médio de dias em um mês. Interessante...

Quaisquer pensamentos ou comentários?

Chris
fonte
0

Pode-se também usar o teste de Ljung-Box para descobrir qual diferença sazonal atinge a melhor estacionariedade. Eu estava trabalhando em um assunto diferente e usei isso realmente para os mesmos fins. Tente períodos diferentes, como 3 a 24, para obter dados mensais. E teste cada um deles pela Ljung-Box e armazene os resultados do Chi-Square. E escolha o período com o menor valor do qui-quadrado.

Aqui está um código simples para fazer isso.

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
todos
fonte