Descomponer los datos de una serie temporal en tendencia determinista y tendencia estocástica

Question

Teniendo en cuenta los comentarios de abajo, he reformado mis preguntas con más información de fondo y simulación/ejemplos. Mi pregunta también incluye la validez del método en teoría. Por favor, corregidme si alguna afirmación no es correcta.

Antecedentes: 1) Para nuestro cliente (los responsables políticos), suele ser interesante comprobar si existe una tendencia lineal en los datos de las series temporales, especialmente tras la aplicación de una nueva política. 2) Una de las características de estos datos es que suelen tener un periodo de tiempo muy corto (unos 10 años). 3) Otra característica de los datos es que la observación no es un censo, sino una estimación a partir de muestras (estudio de biología y ecología de campo, no es posible recoger datos de censo). Esto hace que haya errores de observación en los datos, que a veces pueden ser muy grandes (parecen valores atípicos).

Motivación: Llevo un tiempo intrigado sobre cómo manejar las tres preguntas anteriores. Me gustaría probar todo lo que pueda utilizando ARIMA. Para resolver 1), proporcionando una prueba estadística para una tendencia lineal.

Método: Teóricamente, si los datos son estacionarios después del proceso ARIMA(p,d,q) y contienen una deriva, los datos pueden escribirse como $y_{t}=u\times t + y_{t}^{'}$ y $y_{t}^{'}$ sigue un proceso ARIMA con orden(p,d,q). En un ejemplo de orden (1,1,0), $y_{t}^{'}-y_{t-1}^{'}=\phi(y_{t}^{'}-y_{t-1}^{'})+\epsilon_{t}$ . Estoy pensando en descomponer los datos observados en una tendencia determinista $u\times t $ una tendencia estocástica $y_{t}-u\times t -\epsilon_{t}=y_{t}^{'}-\epsilon_{t}$ y error $\epsilon_{t}$ . ¿Cree que estos dos componentes pueden denominarse en estos términos? Los parámetros de las dos tendencias se pueden estimar mediante la función Arima.

Aquí está el código en un dato simulado para demostrar cómo trabajo a través de él:

set.seed(123)
n  <- 100
e <- rnorm(n,0,1.345)
y1 <- 3.4
AR <- -0.77 
u  <- 0.05
## 1. simulate ARIMA component
ts.sim1 <-arima.sim(n=n,model=list(ar=AR,order=c(1,1,0)),start.innov=y1/(AR),n.start=1,innov=c(0,rnorm(n-1,0,0.345)))
ts.sim1 <- ts(ts.sim1[2:(n+1)])  
ts.sim1
plot(ts.sim1)
## 2. add linear trend
ts.sim2 <- ts.sim1 + u*(1:(n))
plot(ts.sim2)

Este es un trozo de código extra que usé para probar si los parámetros que introduje dan datos estacionarios de ts.sim1 después del proceso (1,1,0).

dat <- replicate(1000, arima.sim(n=n,model=list(ar=AR,order=c(1,1,0)),start.innov=y1/(AR),n.start=1,innov=c(0,rnorm(n-1,0,0.345))))
res  <-  apply(dat, 2, function(x) {fitt <- Arima(x, order=c(1,1,0),  include.drift=F, method="ML"); residuals(fitt)})
p <- apply(res, 2, function(x) adf.test(x)$p.value)
    sum(pv.st > .05)/1000*100

Siguiente:

## 3. make some plots
adf.test(ts.sim2, alternative = "stationary")
Acf(ts.sim2, main='')
Pacf(ts.sim2, main='')

## 4. auto-select best model in terms of AIC, and check residual pattern
fit<-auto.arima(ts.sim2, seasonal=FALSE, trace=TRUE, allowdrift=TRUE)
arima.string1(fit)
tsdisplay(residuals(fit), lag.max=15, main='Best Model Residuals')
AIC(fit)

## 5. Apply a drift version of the best (p,d,q) model, even if the best model does not contain drift, and check residuals
fit1  <-  Arima(ts.sim2, order=c(1,1,0),  include.drift=T, method="ML")
summary(fit1)
AIC(fit);AIC(fit1)  
tsdisplay(residuals(fit1), lag.max=15, main='best Model Residuals')  

Tenga en cuenta que la selección auto.arima no da necesariamente el modelo verdadero, en este caso, sugiere que el mejor modelo es (2,1,0) con deriva en lugar de (1,1,0).

Los gráficos de residuos (abajo) para (1,1,0) con deriva muestran residuos aleatorios independientes con igual varianza

El modelo ARIMA (1,1,0) con deriva arroja el siguiente resultado:

Series: ts.sim2 
ARIMA(1,1,0) with drift         

Coefficients:
          ar1   drift
      -0.9126  0.0298
s.e.   0.0533  0.0192

sigma^2 estimated as 0.1355:  log likelihood=-41.41
AIC=88.83   AICc=89.08   BIC=96.61

Training set error measures:
                       ME      RMSE       MAE       MPE     MAPE      MASE       ACF1
Training set -0.009091779 0.3625225 0.2671725 -3.364276 12.46496 0.5205074 -0.0499297 

A continuación, puedo comprobar la significación estadística de la pendiente lineal

    ## 6. test for linear slope
    drift_index <- 2    
    n <- length(discards)
    pvalue <- 2*pt(-abs(fit1$coef[drift_index]/(sqrt(diag(fit1$var.coef))[drift_index]/sqrt(n))), df=n-1)

Por último, trace la tendencia determinista y estocástica ajustada:

drift_index <- 2
par(mfrow=c(3,1))
t_s <- 1:n
plot(t_s, ts.sim2, type="o", lwd=2, col="red", pch=15, xlab="Year", ylab="", cex.lab=1.5, cex.axis=1.2)

# 1. deterministic trend
ttime <- 1:length(ts.sim2)
y1    <- ttime*fit1$coef[drift_index]
se_re <- sqrt(fit1$sigma2)
m1    <- mean(y1)
offset <- (range(discards)[2]-range(discards)[1])/2 -m1
y1_low <- y1-1.96*se_re+offset
y1_high <- y1+1.96*se_re+offset
plot(t_s,y1+offset, type="n",  ylim=range(y1_low,y1_high), xlab="Year", ylab="Drift", cex.lab=1.5, cex.axis=1.2, main="Deterministic trend", main.cex=1.2)
polygon(c(t_s,rev(t_s)), c(y1_high,rev(y1_low)), 
        col=rgb(0,0,0.6,0.2), border=FALSE)
lines(t_s, y1+offset)

# 2. stochastic trend: fitted value of the ARIMA model part with 95% 
prediction intervals
y2 <- ts.sim2-y1
fitted1 <- y2-residuals(fit1)
se_re <- sqrt(fit1$sigma2)
y2_low <- y2-1.96*se_re
y2_high <- y2+1.96*se_re
plot(t_s, y2, type="n", ylim=range(y2_low,y2_high), xlab="Year", 
ylab="Fitted ARIMA", cex.lab=1.5, cex.axis=1.2, main="Stochastic trend", 
main.cex=1.2)
polygon(c(t_s,rev(t_s)), c(y2_high,rev(y2_low)), 
    col=rgb(0,0,0.6,0.2), border=FALSE)
lines(t_s, y2)

He intentado trazar los intervalos de predicción (condicionales a cada componente) en banda gris en la figura, pero no estoy seguro de haber hecho el cálculo correcto en el código.

Además, he probado una simulación para ver si este proceso de descomposición es insesgado (códigos de abajo). Parece que necesito una serie temporal observada muy larga (n=1000) para obtener una estimación insesgada de los parámetros.

##----simulate N times
N <- 2000
model_type <- rep(NA, N)
res <- data.frame(ar=rep(NA, N), drift=rep(NA,N))
for (k in 1:N) {
  n  <- 1000
  e <- rnorm(n,0,1.345)
  y1 <- 3.4
  AR <- -0.77 
  u  <- 0.05
  ts.sim1 <- arima.sim(n=n,model=list(ar=AR,order=c(1,1,0)),start.innov=y1/(AR),n.start=1,innov=c(0,rnorm(n-1,0,0.345)))
  ts.sim1 <- ts(ts.sim1[2:(n+1)])
  ts.sim2 <- ts.sim1 + u*(1:(n))
  fit<-auto.arima(ts.sim2, seasonal=FALSE, trace=F,allowdrift=TRUE)
  model_type[k] <- arima.string1(fit)
  fit1  <-  Arima(ts.sim2, order=c(1,1,0),  include.drift=T, method="ML")
  res$ar[k] <- coef(fit1)[1]
  res$drift[k] <- coef(fit1)[2]
}
mean(res$ar)
mean(res$drift)
table(model_type)

Mi pregunta es:

1) ¿se está utilizando correctamente la terminología de tendencia determinista frente a la estocástica?

2) Teóricamente, ¿es este un proceso válido para detectar la tendencia lineal y también permitir la observación/error autocorrelacionado? ¿Existe algún otro método para manejar esto bajo ARIMA?

3) En mi última simulación, me he dado cuenta de que la descomposición insesgada sólo funciona cuando la serie temporal es muy larga, lo cual es lo contrario en mis datos (unos 10 años). Supongo que este es el problema general del método ARIMA.

Answer 1

Antes de recibir sus datos, me gustaría tomar el "púlpito" y exponer la tarea que tenemos entre manos y cómo resolvería este enigma. Creo que el enfoque que usted sugiere es formar un modelo ARIMA utilizando procedimientos que implícitamente especifican que no hay variables de tendencia temporal, por lo que se concluye incorrectamente sobre la diferenciación necesaria, etc. Usted asume que no hay valores atípicos, pulsos/impulsos estacionales y que no hay cambios de nivel (cambios de intercepción). Después de una probable especificación errónea del filtro/estructura ARIMA, se asume 1 tendencia/1 intercepción y se reconstruye. Este es un enfoque que, aunque programable, está plagado de defectos lógicos, por no hablar de la varianza de error no constante o de los parámetros no constantes en el tiempo.

El primer paso en el análisis consiste en enumerar el posible espacio muestral que debe investigarse y, en ausencia de una solución directa, realizar una solución basada en el ordenador (ensayo y error) que utiliza una miríada de posibles ensayos/combinaciones que arrojan una posible solución óptima sugerida.

El espacio de la muestra contiene

1. el número de tendencias distintas

2 el número de intercepciones posibles

3 el número y el tipo de operadores de diferenciación

4. la forma del modelo ARMA

5 el número de impulsos únicos

6 el número de pulsos estacionales ( factores estacionales )

7 cualquier punto de cambio de varianza de error requerido, lo que sugiere la necesidad de utilizar mínimos cuadrados ponderados

8 cualquier transformación de potencia necesaria que refleje una relación entre la varianza del error y el valor esperado

Simplemente evalúe todas las permutaciones posibles de estos 8 factores y seleccione aquella combinación única que minimice alguna medida de error porque ¡el orden es importante! .

Si esto es oneroso, que así sea y espero recibir su tsim2 para poder (posiblemente) demostrar un enfoque que hable de este "espinoso tema" usando algunos de mis juguetes favoritos.

Tenga en cuenta que si ha simulado (estrechamente) entonces su enfoque podría ser la respuesta, pero la pregunta que tengo es "su enfoque es robusto a las violaciones de datos" o es simplemente un enfoque de libro de cocina que funciona en este conjunto de datos y falla en otros. Confía pero verifica.

EDITADO TRAS LA RECEPCIÓN DE LOS DATOS (100 VALORES)

Confío en que esta discusión ponga de manifiesto la necesidad de enfoques exhaustivos/programables para formar modelos útiles. Como se ha comentado anteriormente, un torneo eficiente basado en el ordenador, en el que se examinaron las diferentes combinaciones posibles (un máximo de 256), dio como resultado la siguiente propuesta de modelo inicial.

El concepto aquí es "duplicar/aproximar el ojo humano" examinando las alternativas que compiten, que es lo que (en mi opinión) hacemos cuando realizamos la identificación visual de la estructura. En este caso, la mayoría de los ojos no verán el cambio de nivel en el período 65 y simplemente se centrarán en la principal ruptura de la tendencia en torno al período 51.

1 IDENTIFICAR PUNTOS DE RUPTURA DETERMINISTAS EN LA TENDENCIA 2 IDENTIFICAR LOS CAMBIOS DE INTERCEPCIÓN 2 EVALUAR LA NECESIDAD DE AUMENTAR LA ARIMA 4 EVALUAR LA NECESIDAD DE PULSOS

SIMPLIFY VIA NECESSITY TESTS 

detallando tanto un cambio de tendencia (51) como un cambio de intercepción (65). La comprobación del diagnóstico del modelo (siempre una buena idea en los enfoques iterativos de la forma del modelo) dio como resultado la siguiente acf sugiriendo que era necesaria una mejora para obtener un conjunto de residuos sin estructura. Se propuso entonces un modelo aumentado de la forma con un coeficiente AR(1) insignificante.

El modelo final está aquí con las estadísticas del modelo y aquí

Los residuos de este modelo se presentan aquí con un acf de

El gráfico de la realidad/ajuste y la previsión está aquí . El limpiado frente al actual es revelador ya que detalla el efecto de cambio de nivel

En resumen, donde el OP simuló un (1,1,0) para las primeras 50 observaciones, luego abrevió las últimas 50 observaciones coloreando/cambiando efectivamente el proceso ARMA compuesto a un (1,0,0) mientras incorpora los 3 predictores identificados empíricamente.

El objetivo es un análisis exhaustivo de los datos que incorpore procedimientos de búsqueda avanzados. Este conjunto de datos es "espinoso" y espero cualquier sugerencia de mejora que pueda surgir de este debate. He utilizado una versión beta de AUTOBOX (que he ayudado a desarrollar) como herramienta de elección.

En cuanto a su "método propuesto", puede funcionar para esta serie, pero hay demasiadas suposiciones, como una y sólo una tendencia estocástica, una y sólo una tendencia determinista (1,2,3,...), sin pulsos, sin cambios de nivel (cambios de intercepción), sin pulsos estacionales, varianza de error constante, parámetros constantes a lo largo del tiempo, etc., para sugerir la generalidad del enfoque. Usted está argumentando de lo específico a lo general. Hay toneladas de soluciones ad hoc erróneas esperando a ser especificadas y sólo un puñado de "soluciones correctas" de las cuales mi enfoque es sólo una.

Un primer plano que muestra las observaciones 51 a 100 sugiere una desviación/cambio significativo en el patrón (es decir, intercepción implícita) a partir del periodo 65 ( que fue recogido/identificado por el análisis como un cambio de nivel (cambio de intercepción)) lo que sugiere un posible fallo de simulación, ya que las obs. 51-64 tienen un patrón diferente al de las obs. 65-100.