tendencia/estacionalidad estocástica frente a determinista en la previsión de series temporales

Question

Tengo una formación moderada en previsión de series temporales. He consultado varios libros de previsión y no veo que en ninguno de ellos se aborden las siguientes cuestiones.

Tengo dos preguntas:

1. ¿Cómo podría determinar objetivamente (mediante una prueba estadística) si una serie temporal dada tiene:

   • Estacionalidad estocástica o estacionalidad determinista
   • Tendencia estocástica o tendencia determinista
2. ¿Qué pasaría si modelo mi serie temporal como una tendencia/estacionalidad determinista cuando la serie tiene un componente claramente estocástico?

Cualquier ayuda para responder a estas preguntas será muy apreciada.

Ejemplo de datos para la tendencia:

7,657
5,451
10,883
9,554
9,519
10,047
10,663
10,864
11,447
12,710
15,169
16,205
14,507
15,400
16,800
19,000
20,198
18,573
19,375
21,032
23,250
25,219
28,549
29,759
28,262
28,506
33,885
34,776
35,347
34,628
33,043
30,214
31,013
31,496
34,115
33,433
34,198
35,863
37,789
34,561
36,434
34,371
33,307
33,295
36,514
36,593
38,311
42,773
45,000
46,000
42,000
47,000
47,500
48,000
48,500
47,000
48,900

Answer 1

1) En cuanto a su primera pregunta, en la literatura se han desarrollado y discutido algunas pruebas estadísticas para probar la nulidad de estacionariedad y la nulidad de raíz unitaria. Algunos de los muchos trabajos que se han escrito sobre este tema son los siguientes:

Relacionado con la tendencia:

• Dickey, D. y Fuller, W. (1979a), Distribución de los estimadores de autoregressive time series with a unit root, Journal of the American Statistical Association 74, 427-31.
• Dickey, D. y Fuller, W. (1981), Likelihood ratio statistics for autoregressive time series with a unit root, Econometrica 49, 1057-1071.
• Kwiatkowski, D., Phillips, P., Schmidt, P. y Shin, Y. (1992), Testing the null hypothesis of stationarity against the alternative of a unit root: How sure Journal of Econometrics, 54, 159-178. 54, 159-178.
• Phillips, P. y Perron, P. (1988), Testing for a unit root in time series regression, Biometrika 75, 335-46.
• Durlauf, S. y Phillips, P. (1988), Trends versus random walks in time series analysis, Econometrica 56, 1333-54.

Relacionado con el componente estacional:

• Hylleberg, S., Engle, R., Granger, C. y Yoo, B. (1990), Integración estacional and cointegration, Journal of Econometrics 44, 215-38.
• Canova, F. y Hansen, B. E. (1995), ¿Son constantes en el tiempo los patrones estacionales? a test for seasonal stability, Journal of Business and Economic Statistics 13, 237-252.
• Franses, P. (1990), Testing for seasonal unit roots in monthly data, Informe técnico 9032, Econometric Institute.
• Ghysels, E., Lee, H. y Noh, J. (1994), Testing for unit roots in seasonal time series. some theoretical extensions and a monte carlo investigation, Journal of Econometrics 62, 415-442.

El libro de texto Banerjee, A., Dolado, J., Galbraith, J. y Hendry, D. (1993), Co-Integration,Error Correction, and the econometric analysis of non-stationary data, Advanced Texts in Econometrics. Oxford University Press es también una buena referencia.

2) Su segunda preocupación está justificada por la bibliografía. Si existe una prueba de raíz unitaria, el estadístico t tradicional que se aplicaría a una tendencia lineal no sigue la distribución estándar. Véase, por ejemplo, Phillips, P. (1987), Time series regression with unit root, Econometrica 55(2), 277-301.

Si existe una raíz unitaria y se ignora, entonces se reduce la probabilidad de rechazar la nulidad de que el coeficiente de una tendencia lineal sea cero. Es decir, acabaríamos modelizando una tendencia lineal determinista con demasiada frecuencia para un nivel de significación dado. En presencia de una raíz unitaria, en su lugar deberíamos transformar los datos tomando diferencias regulares a los datos.

3) A modo ilustrativo, si utilizas R puedes hacer el siguiente análisis con tus datos.

x <- structure(c(7657, 5451, 10883, 9554, 9519, 10047, 10663, 10864, 
  11447, 12710, 15169, 16205, 14507, 15400, 16800, 19000, 20198, 
  18573, 19375, 21032, 23250, 25219, 28549, 29759, 28262, 28506, 
  33885, 34776, 35347, 34628, 33043, 30214, 31013, 31496, 34115, 
  33433, 34198, 35863, 37789, 34561, 36434, 34371, 33307, 33295, 
  36514, 36593, 38311, 42773, 45000, 46000, 42000, 47000, 47500, 
  48000, 48500, 47000, 48900), .Tsp = c(1, 57, 1), class = "ts")

En primer lugar, puede aplicar la prueba Dickey-Fuller para la nulidad de una raíz unitaria:

require(tseries)
adf.test(x, alternative = "explosive")
#   Augmented Dickey-Fuller Test
#   Dickey-Fuller = -2.0685, Lag order = 3, p-value = 0.453
#   alternative hypothesis: explosive

y la prueba KPSS para la hipótesis nula inversa, estacionariedad frente a la alternativa de estacionariedad en torno a una tendencia lineal:

kpss.test(x, null = "Trend", lshort = TRUE)
#   KPSS Test for Trend Stationarity
#   KPSS Trend = 0.2691, Truncation lag parameter = 1, p-value = 0.01

Resultados: Prueba ADF, al nivel de significación del 5% no se rechaza una raíz unitaria; Prueba KPSS, se rechaza el nulo de estacionariedad a favor de un modelo con tendencia lineal.

Nota al margen: el uso de lshort=FALSE no se rechaza la nulidad de la prueba KPSS al nivel del 5%, sin embargo, selecciona 5 rezagos; una inspección posterior no mostrada aquí sugiere que la elección de 1-3 rezagos es apropiada para los datos y lleva a rechazar la hipótesis nula.

En principio, deberíamos guiarnos por la prueba para la que pudimos rechazar la hipótesis nula (en lugar de por la prueba para la que no rechazamos (aceptamos) la nula). Sin embargo, una regresión de la serie original sobre una tendencia lineal resulta no ser fiable. Por un lado, la R-cuadrado es elevado (más del 90%), lo que se señala en la bibliografía como un indicador de regresión espuria.

fit <- lm(x ~ 1 + poly(c(time(x))))
summary(fit)
#Coefficients:
#                 Estimate Std. Error t value Pr(>|t|)    
#(Intercept)       28499.3      381.6   74.69   <2e-16 ***
#poly(c(time(x)))  91387.5     2880.9   31.72   <2e-16 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#Residual standard error: 2881 on 55 degrees of freedom
#Multiple R-squared:  0.9482,   Adjusted R-squared:  0.9472 
#F-statistic:  1006 on 1 and 55 DF,  p-value: < 2.2e-16

Por otra parte, los residuos están autocorrelacionados:

acf(residuals(fit)) # not displayed to save space

Además, no se puede rechazar la nulidad de una raíz unitaria en los residuos.

adf.test(residuals(fit))
#   Augmented Dickey-Fuller Test
#Dickey-Fuller = -2.0685, Lag order = 3, p-value = 0.547
#alternative hypothesis: stationary

En este punto, se puede elegir el modelo que se utilizará para obtener las previsiones. Por ejemplo, las previsiones basadas en un modelo estructural de series temporales y en un modelo ARIMA pueden obtenerse del siguiente modo.

# StructTS
fit1 <- StructTS(x, type = "trend")
fit1
#Variances:
# level    slope  epsilon  
#2982955        0   487180 
# 
# forecasts
p1 <- predict(fit1, 10, main = "Local trend model")
p1$pred
# [1] 49466.53 50150.56 50834.59 51518.62 52202.65 52886.68 53570.70 54254.73
# [9] 54938.76 55622.79

# ARIMA
require(forecast)
fit2 <- auto.arima(x, ic="bic", allowdrift = TRUE)
fit2
#ARIMA(0,1,0) with drift         
#Coefficients:
#         drift
#      736.4821
#s.e.  267.0055
#sigma^2 estimated as 3992341:  log likelihood=-495.54
#AIC=995.09   AICc=995.31   BIC=999.14
#
# forecasts
p2 <- forecast(fit2, 10, main = "ARIMA model")
p2$mean
# [1] 49636.48 50372.96 51109.45 51845.93 52582.41 53318.89 54055.37 54791.86
# [9] 55528.34 56264.82

Un gráfico de las previsiones:

par(mfrow = c(2, 1), mar = c(2.5,2.2,2,2))
plot((cbind(x, p1$pred)), plot.type = "single", type = "n", 
  ylim = range(c(x, p1$pred + 1.96 * p1$se)), main = "Local trend model")
grid()
lines(x)
lines(p1$pred, col = "blue")
lines(p1$pred + 1.96 * p1$se, col = "red", lty = 2)
lines(p1$pred - 1.96 * p1$se, col = "red", lty = 2)
legend("topleft", legend = c("forecasts", "95% confidence interval"), 
  lty = c(1,2), col = c("blue", "red"), bty = "n")
plot((cbind(x, p2$mean)), plot.type = "single", type = "n", 
  ylim = range(c(x, p2$upper)), main = "ARIMA (0,1,0) with drift")
grid()
lines(x)
lines(p2$mean, col = "blue")
lines(ts(p2$lower[,2], start = end(x)[1] + 1), col = "red", lty = 2)
lines(ts(p2$upper[,2], start = end(x)[1] + 1), col = "red", lty = 2)

Las previsiones son similares en ambos casos y parecen razonables. Observe que las previsiones siguen un patrón relativamente determinista similar a una tendencia lineal, pero no modelizamos explícitamente una tendencia lineal. La razón es la siguiente i) en el modelo de tendencia local, la varianza del componente de pendiente se estima como cero. Esto convierte el componente de tendencia en una deriva que tiene el efecto de una tendencia lineal. ii) ARIMA(0,1,1), se selecciona un modelo con una deriva en un modelo para las series diferenciadas.El efecto del término constante en una serie diferenciada es una tendencia lineal. Este tema se trata en esta entrada .

Puede comprobar que si se elige un modelo local o un ARIMA(0,1,0) sin deriva, entonces las previsiones son una línea recta horizontal y, por tanto, no tendrían ningún parecido con la dinámica observada de los datos. Pues bien, esto forma parte del rompecabezas de las pruebas de raíz unitaria y los componentes deterministas.

Edición 1 (inspección de los residuos): La autocorrelación y la ACF parcial no sugieren una estructura en los residuos.

resid1 <- residuals(fit1)
resid2 <- residuals(fit2)
par(mfrow = c(2, 2))
acf(resid1, lag.max = 20, main = "ACF residuals. Local trend model")
pacf(resid1, lag.max = 20, main = "PACF residuals. Local trend model")
acf(resid2, lag.max = 20, main = "ACF residuals. ARIMA(0,1,0) with drift")
pacf(resid2, lag.max = 20, main = "PACF residuals. ARIMA(0,1,0) with drift")

Como sugirió IrishStat, también es aconsejable comprobar la presencia de valores atípicos. Se detectan dos valores atípicos aditivos utilizando el paquete tsoutliers .

require(tsoutliers)
resol <- tsoutliers(x, types = c("AO", "LS", "TC"), 
  remove.method = "bottom-up", 
  args.tsmethod = list(ic="bic", allowdrift=TRUE))
resol
#ARIMA(0,1,0) with drift         
#Coefficients:
#         drift        AO2       AO51
#      736.4821  -3819.000  -4500.000
#s.e.  220.6171   1167.396   1167.397
#sigma^2 estimated as 2725622:  log likelihood=-485.05
#AIC=978.09   AICc=978.88   BIC=986.2
#Outliers:
#  type ind time coefhat  tstat
#1   AO   2    2   -3819 -3.271
#2   AO  51   51   -4500 -3.855

Observando la ACF, podemos decir que, al nivel de significación del 5%, los residuos son aleatorios también en este modelo.

par(mfrow = c(2, 1))
acf(residuals(resol$fit), lag.max = 20, main = "ACF residuals. ARIMA with additive outliers")
pacf(residuals(resol$fit), lag.max = 20, main = "PACF residuals. ARIMA with additive outliers")

En este caso, la presencia de posibles valores atípicos no parece distorsionar el rendimiento de los modelos. Así lo corrobora la prueba de normalidad de Jarque-Bera; el nulo de normalidad en los residuos de los modelos iniciales ( fit1 , fit2 ) no se rechaza al nivel de significación del 5%.

jarque.bera.test(resid1)[[1]]
# X-squared = 0.3221, df = 2, p-value = 0.8513
jarque.bera.test(resid2)[[1]]
#X-squared = 0.426, df = 2, p-value = 0.8082

Edición 2 (gráfico de los residuos y sus valores) Este es el aspecto de los residuos:

Y estos son sus valores en formato csv:

0;6.9205
-0.9571;-2942.4821
2.6108;4695.5179
-0.5453;-2065.4821
-0.2026;-771.4821
0.1242;-208.4821
0.1909;-120.4821
-0.0179;-535.4821
0.1449;-153.4821
0.484;526.5179
1.0748;1722.5179
0.3818;299.5179
-1.061;-2434.4821
0.0996;156.5179
0.4805;663.5179
0.8969;1463.5179
0.4111;461.5179
-1.0595;-2361.4821
0.0098;65.5179
0.5605;920.5179
0.8835;1481.5179
0.7669;1232.5179
1.4024;2593.5179
0.3785;473.5179
-1.1032;-2233.4821
-0.3813;-492.4821
2.2745;4642.5179
0.2935;154.5179
-0.1138;-165.4821
-0.8035;-1455.4821
-1.2982;-2321.4821
-1.9463;-3565.4821
-0.1648;62.5179
-0.1022;-253.4821
0.9755;1882.5179
-0.5662;-1418.4821
-0.0176;28.5179
0.5;928.5179
0.6831;1189.5179
-1.8889;-3964.4821
0.3896;1136.5179
-1.3113;-2799.4821
-0.9934;-1800.4821
-0.4085;-748.4821
1.2902;2482.5179
-0.0996;-657.4821
0.5539;981.5179
2.0007;3725.5179
1.0227;1490.5179
0.27;263.5179
-2.336;-4736.4821
1.8994;4263.5179
0.1301;-236.4821
-0.0892;-236.4821
-0.1148;-236.4821
-1.1207;-2236.4821
0.4801;1163.5179

Answer 2

Con respecto a sus datos no estacionales ...Las tendencias pueden ser de dos formas y(t)=y(t-1)+θ0 (A) Tendencia estocástica o Y(t)=a+bx1+cx2 (B) Tendencia determinista etc donde x1=1,2,3,4....t y x2=0,0,0,0,0,1,2,3,4 por lo que una tendencia se aplica a las observaciones 1-t y una segunda tendencia se aplica a las observaciones 6 a t.

Su serie no estacional contenía 29 valores. Utilicé AUTOBOX, un programa informático que había ayudado a desarrollar de forma totalmente automática. AUTOBOX es un procedimiento transparente, ya que detalla cada paso del proceso de modelización. Aquí se presenta un gráfico de las series/valores ajustados/previsiones . La utilización de AUTOBOX para formar un modelo de tipo A condujo a lo siguiente . La ecuación se presenta de nuevo aquí Los estadísticos del modelo son . Aquí se muestra un gráfico de los residuos mientras que la tabla de valores previstos está aquí . Al restringir AUTOBOX a un modelo de tipo B, AUTOBOX detectó una tendencia creciente en el periodo 14:. ¡!

En cuanto a la comparación de modelos: Dado que el número de observaciones ajustadas difiere (26 y 29 respectivamente), no es posible utilizar métricas estándar (es decir, r-cuadrado, desviación estándar del error, AIC, etc.) para determinar la dominancia, aunque en este caso el asentimiento sería para A. Los residuos de A son mejores debido a la estructura AR(2). Las previsiones de B son un poco agresivas, mientras que el patrón de las previsiones de A es más intuitivo. Se podrían retener, por ejemplo, 4 observaciones y evaluar la precisión de una previsión a 1 periodo vista a partir de 4 orígenes distintos (25, 26, 27 y 28).