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

Esta publicación 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 punto estoy interesado en una forma sólida de descubrir la periodicidad / estacionalidad de una serie temporal genérica afectada por mucho ruido. Desde el punto de vista del desarrollador, me gustaría una interfaz simple como:

unsigned int discover_period(vector<double> v);

Dónde vestá la matriz que contiene las muestras y el valor de retorno es el período de la señal. El punto principal es que, nuevamente, no puedo hacer ninguna suposición con respecto a la señal analizada. Ya probé un enfoque basado en la autocorrelación de señal (detectando los picos de un correlograma), pero no es robusto como me gustaría.

Si realmente no tiene idea de cuál es la periodicidad, probablemente el mejor enfoque es encontrar la frecuencia correspondiente al máximo de la densidad espectral. Sin embargo, el espectro a bajas frecuencias se verá afectado por la tendencia, por lo que primero debe reducir la tendencia de la serie. La siguiente función R debería hacer el trabajo para la mayoría de las series. Está lejos de ser perfecto, pero lo he probado en algunas docenas de ejemplos y parece funcionar bien. Devolverá 1 para datos que no tienen una periodicidad fuerte, y la duración del período de lo contrario.

Actualización: 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)
}

Si espera que el proceso sea estacionario (la periodicidad / estacionalidad no cambiará con el tiempo), entonces algo como un periodograma Chi-cuadrado (ver, por ejemplo, Sokolove y Bushell, 1978) podría ser una buena opción. Se usa comúnmente en el análisis de datos circadianos que pueden tener cantidades extremadamente grandes de ruido, pero se espera que tenga periodicidades muy estables.

Este enfoque no presupone la forma de la forma de onda (aparte de que es consistente de ciclo a ciclo), pero requiere que cualquier ruido sea de media constante y no esté correlacionado con 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 solo un ejemplo, y muestran que puede identificar el período de una función trigonométrica pura, incluso con mucho ruido aditivo.

Tal como está escrito, el último argumento ( alpha) en la llamada es superfluo, la función simplemente devuelve el "mejor" período que puede encontrar; descomente la primera returndeclaración y comente la segunda para que devuelva una lista de todos los períodos significativos en el nivel alpha.

Esta función no realiza ningún tipo de comprobación de la cordura para asegurarse de que ha puesto en períodos identificables, ni funciona (puede) con períodos fraccionarios, ni hay ningún tipo de control de comparación múltiple incorporado si decide mira múltiples períodos. Pero aparte de eso, debería ser razonablemente robusto.

Es posible que desee definir lo que quiere más claramente (para usted mismo, si no aquí). Si lo que está buscando es el período estacionario estadísticamente más significativo contenido en sus datos ruidosos, esencialmente hay dos rutas a seguir:

1) calcule una estimación de autocorrelación robusta y tome el coeficiente máximo
2) calcule una estimación de densidad espectral de potencia robusta y tome el máximo del espectro

El problema con el n. ° 2 es que para cualquier serie temporal ruidosa, obtendrá una gran cantidad de potencia en bajas frecuencias, lo que hace que sea difícil distinguirla. Existen algunas técnicas para resolver este problema (es decir, pre-blanqueamiento, luego estimar el PSD), pero si el período verdadero de sus datos es lo suficientemente largo, la detección automática será dudosa.

Su mejor opción es probablemente implementar una rutina de autocorrelación robusta como se puede encontrar en el capítulo 8.6, 8.7 en Estadísticas robustas - Teoría y métodos de Maronna, Martin y Yohai. La búsqueda en Google de "robusto durbin-levinson" también arrojará algunos resultados.

Si solo está buscando una respuesta simple, no estoy seguro de que exista. La detección de períodos en series de tiempo puede ser complicada, y pedir una rutina automatizada que pueda realizar magia puede ser demasiado.

Podría usar la transformación de Hilbert de la teoría DSP para medir la frecuencia instantánea de sus datos. El sitio http://ta-lib.org/ tiene un código fuente abierto para medir el período del ciclo dominante de datos financieros; la función relevante se llama HT_DCPERIOD; es posible que pueda usar esto o adaptar el código a sus propósitos.

Un enfoque diferente podría ser la descomposición en modo empírico. El paquete R se llama EMD desarrollado por el inventor del 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()

El método fue calificado como 'Empírico' por una buena razón y existe el riesgo de que las Funciones del Modo Intrínseco (los componentes aditivos individuales) se mezclen. Por otro lado, el método es muy intuitivo y puede ser útil para una inspección visual rápida de la ciclicidad.

En referencia a la publicación de Rob Hyndman arriba https://stats.stackexchange.com/a/1214/70282

La función find.freq funciona de manera brillante. En el conjunto de datos diarios que estoy usando, funcionó correctamente la frecuencia para ser 7.

Cuando lo probé solo los días de la semana, mencionó que la frecuencia es 23, que es notablemente cercana a 21.42857 = 29.6 * 5/7, que es el número promedio de días de trabajo en un mes. (O, por el contrario, 23 * 7/5 es 32).

Mirando hacia atrás a mis datos diarios, experimenté con el presentimiento de tomar el primer período, promediar con eso y luego encontrar el siguiente período, etc. Ver a continuación:

find.freq.all = function (x) {  
  f = find.freq (x);
  freqs = c (f);  
  mientras que (f> 1) {
    inicio = 1; #también intente iniciar = f;
    x = period.apply (x, seq (inicio, longitud (x), f), media); 
    f = find.freq (x);
    freqs = c (freqs, f);
  }
  if (length (freqs) == 1) {return (freqs); }
  para (i en 2: longitud (freqs)) {
    freqs [i] = freqs [i] * freqs [i-1];
  }
  freqs [1: (longitud (freqs) -1)];
}
find.freq.all (dailyts) #usando datos diarios

Lo anterior da (7,28) o (7,35) dependiendo de si la secuencia comienza con 1 o f. (Ver comentario arriba)

Lo que implicaría que los períodos estacionales para msts (...) deberían 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 es 31.5, que está cerca de la duración promedio de un mes.

Sospecho que reinventé la rueda, ¿cómo se llama este algoritmo? ¿Hay una mejor implementación en R en alguna parte?

Más tarde, ejecuté el código anterior al intentar todos los inicios del 1 al 7 y obtuve 35,35,28,28,28,28,28 para el segundo período. El promedio es de 30, que es el número promedio de días en un mes. Interesante...

¿Alguna idea o comentario?

También se puede usar la prueba de Ljung-Box para descubrir qué diferencia estacional alcanza la mejor estacionariedad. Estaba trabajando en un tema diferente y en realidad lo usé para los mismos propósitos. Pruebe diferentes períodos, como 3 a 24 para obtener datos mensuales. Y pruebe cada uno de ellos con Ljung-Box y almacene los resultados de Chi-Square. Y elija el período con el menor valor de chi-cuadrado.

Aquí hay un código simple para hacer eso.

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