Período detección de una serie temporal genérico

Question

Este post es la continuación de otra publicación relacionada con un método genérico para la detección de valores atípicos en series de tiempo. Básicamente, en este momento estoy interesado en un camino sólido para descubrir la periodicidad/estacionalidad de la genérica de la serie de tiempo de los afectados por una gran cantidad de ruido. De un desarrollador punto de vista, me gustaría una simple interfaz, tales como:

unsigned int discover_period(vector v);

Donde v es la matriz que contiene las muestras, y el valor de retorno es el periodo de la señal. El punto principal es que, de nuevo, no puedo hacer ninguna suposición acerca de la señal analizada. Ya lo he intentado un enfoque basado en la señal de autocorrelación (la detección de los picos de un correlogram), pero no es robusta como me gustaría.

Answer 1

Si realmente no tienes idea de lo que la periodicidad es, probablemente, el mejor enfoque es encontrar la frecuencia que corresponde al máximo de la densidad espectral. Sin embargo, el espectro en frecuencias bajas se verán afectadas por la tendencia, por lo que necesita para detrend la primera serie. La siguiente función de R debe hacer el trabajo para la mayoría de las series. Es lejos de ser perfecto, pero lo he probado en un par de docenas de ejemplos y parece que funciona ok. Se devolverá 1 para los datos que no tienen una fuerte periodicidad y la duración del período de otra manera.

Actualización: la Versión 2 de la función. Esto es mucho más rápido y parece ser más 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)
}

Answer 2

Si usted espera que el proceso sea estacionario -- la periodicidad/estacionalidad no va a cambiar a lo largo del tiempo, a continuación, algo así como una Chi-cuadrado periodograma (ver, por ejemplo, Sokolove y Bushell, 1978) podría ser una buena opción. Es comúnmente usado en el análisis de circadiano de datos que pueden tener cantidades muy grandes de ruido, pero se espera que muy estable periodicidades.

Este enfoque no hace ninguna suposición acerca de la forma de la onda (que es consistente de ciclo a ciclo), pero no requiere que cualquier ruido constante de la media y la correlación de la señal.

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)

Las dos últimas líneas son sólo un ejemplo, mostrando que se puede identificar el período puro de una función trigonométrica, incluso con un montón de ruido aditivo.

Como está escrito, el último argumento (alfa) en la llamada es superfluo, la función simplemente devuelve el 'mejor' período se puede encontrar; elimine la primera de "devolución" de la declaración y comentar el segundo a tener que devolver una lista de todos los períodos significativos en el nivel alfa.

Esta función no hace ningún tipo de cordura comprobación para asegurarse de que usted ha puesto en períodos de identificación personal, ni tampoco (puede) el trabajo con fracciones de los períodos, ni hay ningún tipo de comparación múltiple integrado de control de si usted decide buscar en varios períodos. Pero aparte de eso debería ser razonablemente robusto.

Answer 3

Puede utilizar la transformación de Hilbert de la teoría DSP para medir la frecuencia instantánea de los datos. El sitio http://ta-lib.org/ tiene código fuente abierto para medir el periodo dominante del ciclo de información financiera; la función correspondiente se llama HT_DCPERIOD; usted podría ser capaz de usar este o adaptar el código a sus propósitos.

Answer 4

Es posible que desee definir lo que desea más claramente (a ti mismo, si no aquí). Si lo que estás buscando es el más significativo estadísticamente estacionaria período contenida en sus datos ruidosos, hay básicamente dos caminos a tomar:

1) calcular una robusta de autocorrelación de la estimación, y aprovechar al máximo el coeficiente de
2) calcular un sólido de densidad espectral de potencia de la estimación, y tomar el máximo del espectro

El problema con el #2 es que para cualquier ruidoso de series de tiempo, usted recibirá una gran cantidad de energía en las bajas frecuencias, lo que hace difícil distinguir. Hay algunas técnicas para la resolución de este problema (es decir, pre-blanquear, entonces la estimación de la PSD), pero si el verdadero período de sus datos es lo suficientemente largo, la detección automática será dudoso.

Su mejor apuesta es, probablemente, la implementación de un completo autocorrelación de rutina, tales como pueden ser encontrados en el capítulo 8.6, 8.7 en Estadísticas Robustas - Teoría y Métodos por Maronna, Martin y Yohai. Buscar en Google por "robusto de durbin-levinson" también a producir algunos resultados.

Si usted está buscando una respuesta simple, no estoy seguro de que uno existe. Período de detección en las series de tiempo puede ser complicado, y pidiendo una rutina automatizada que puede realizar magia puede ser demasiado.

Answer 5

En referencia a Rob Hyndman el post de arriba http://stats.stackexchange.com/a/1214/70282

El hallazgo.freq función funciona a la perfección. En el conjunto de datos que estoy utilizando, correctamente trabajado la frecuencia de 7.

Cuando me lo probé sólo los días de la semana, se menciona la frecuencia es de 23, que está muy cerca de 21.42857=29.6*5/7, que es el número promedio de días de trabajo en un mes. (O a la inversa, 23*7/5 es de 32.)

Mirando hacia atrás en mis datos diarios, experimenté con una corazonada de tomar el primer período, con un promedio de por que y, a continuación, encontrar el próximo período, etc. Ver a continuación:

encontrar.freq.todos=function(x){ 
f=buscar.freq(x);
 freqs=c(f); 
mientras(f>1){
 inicio=1; #también intentar start=f;
 x=período.aplica(x,seq(inicio,longitud(x),f),media); 
f=buscar.freq(x);
freqs=c(freqs,f);
}
 si(largo(freqs)==1){ return(freqs); }
 for(i in 2:longitud(freqs)){
freqs[i]=freqs[i]*freqs[i-1];
}
freqs[1:(length(freqs)-1)];
}
encontrar.freq.todos(dailyts) #utilizando datos diarios

El de arriba te da (7,28) o (7,35) dependiendo de si la secuencia se inicia con 1 o f. (Ver comentario anterior).

Lo cual implicaría que los períodos estacionales para ha(...) debe ser (7,28) o (7,35).

La lógica parece sensible a las condiciones iniciales, dada la sensibilidad de los parámetros del algoritmo. La media de 28 y 35 años es del 31,5 que está cerca de la longitud media de un mes.

Sospecho que me reinventado la rueda, ¿cuál es el nombre de este algoritmo? Existe una mejor implementación en R en algún lugar?

Más tarde, me encontré con el código de arriba en tratar todo comienza de 1 a 7 y tengo 35,35,28,28,28,28,28 para el segundo periodo. El promedio de las obras a 30, que es el número promedio de días en un mes. Interesante...

Cualquier pensamiento o comentarios?