Tengo datos de series temporales que necesito transformar en una serie estacionaria. Se trata de una serie mensual, por lo que, según este gráfico, creo que el punto de ruptura está en septiembre de 2007. Sin embargo, no estoy seguro de si debo seguir mi corazonada o encontrar un módulo en R para obtener el punto de ruptura exacto.
Respuestas
¿Demasiados anuncios?Hay muchos paquetes de R. He aquí un resumen no exhaustivo con ejemplos prácticos . A continuación, mostraré una solución utilizando la función mcp
paquete.
Empecemos por simular unos datos que se parezcan un poco a los tuyos:
df = data.frame(x = 1:150)
df$y = c(arima.sim(n = 75, list(ar = 0.8), sd = 2) + 0,
arima.sim(n = 75, list(ar = 0.5), sd = 2) + 7)
Puede ajustar un modelo consistente en dos mesetas con residuos autorregresivos de primer orden.
model = list(
y ~ 1 + ar(1),
~ 1 + ar(1)
)
Tal vez sólo desee un coeficiente autorregresivo compartido para ambos segmentos. En ese caso, elimine ar(1)
en la fórmula para el segundo segmento - que hará que el primer ar(1)
traspasar.
Ahora ajustamos el modelo:
library(mcp)
fit = mcp(model, data = df, par_x = "x")
A veces quieres hacer más muestreos y más rápido a través de mcp(..., iter = 10000, cores = 3)
.
Veamos el resultado visualmente:
plot(fit)
Parece que captura bien los datos. Los puntos negros son los datos brutos. Las líneas grises son muestras de la distribución posterior. La distribución azul es el punto de cambio inferido. En este caso, el cambio en la intercepción es probablemente la "señal" más clara, pero este artículo sobre análisis de series temporales con mcp muestra ejemplos sin cambios tan bruscos en la media.
Puede ver los valores posteriores de cada parámetro utilizando plot_pars(fit)
. Veamos también si recuperó bien los parámetros:
> summary(fit)
Population-level parameters:
name mean lower upper Rhat n.eff
ar1_1 0.74 0.56 0.91 1.0 7131
ar1_2 0.61 0.43 0.80 1.0 6133
cp_1 75.48 72.05 76.00 1.2 147
int_1 0.64 -1.11 2.50 1.0 6288
int_2 6.56 5.22 7.97 1.0 4948
sigma_1 2.09 1.86 2.35 1.0 6743
Simulamos los datos de ar1_1 = 0.8
, ar1_2 = 0.5
, cp_1 = 75
, int_1 = 0
, int_2 = 7
y sigma_1 = 2
. Así que se ve bien.
Se trataba de un modelo muy simple con sólo interceptos y un único punto de cambio. Varios otros paquetes son adecuados aquí también, incluyendo segmented
y EnvCpt
. Puntos fuertes de mcp
incluyen (1) cuantifica la incertidumbre de los puntos de cambio, (2) puede hacer modelos de regresión por segmento, y (3) puede manejar estructuras autorregresivas más complejas, entre otras cosas. Las principales limitaciones son (1) es más lento, pero no excesivamente, (2) requiere JAGS (que es gratuito).
El paquete "changepoint package in R" antes mencionado supone que no hay memoria arima, es decir, que los datos son independientes y se distribuyen normalmente véase https://www.jstatsoft.org/article/view/v058i03 sección 4.1 . Si los datos no son independientes, habrá que identificar el modelo arima subyacente, pero no se ha previsto la identificación de este componente crítico, sino que sólo se pide al usuario que lo especifique previamente.
Además, si el usuario tiene alguna noción sobre los predictores causales y, posiblemente, sus efectos contemporáneos y retardados que podrían necesitar ser condicionados para no mucho lo que limita su utilidad general.
Los breakpoints pueden producirse en el componente arima a lo largo del tiempo y, ciertamente, también podría ser necesario condicionar los puntos de cambio de la varianza de error para ver http://docplayer.net/12080848-Outliers-level-shifts-and-variance-changes-in-time-series.html . Hay que CONOCER los supuestos en los que son válidas todas y cada una de las soluciones propuestas. No es oro todo lo que reluce.