Este proyecto trata sobre ajustar un modelo a la serie de tiempo del precio promedio mensual del azúcar en la Central de Abastos del Distrito Federal. La serie es del precio del azúcar stándar por bulto (50kg) y empieza en octubre del 2000 y termina en noviembre del 2013.
Cargamos las bibliotecas para el análisis de series de tiempo en R.
library(reshape)
## Loading required package: plyr
##
## Attaching package: 'reshape'
##
## The following objects are masked from 'package:plyr':
##
## rename, round_any
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(tseries)
library(knitr)
Primero leemos la base de datos, obtenemos el promedio por mes-mercado-tipo y escogemos el precio de la Central de Abastos del DF para azúcar stándar. Posteriormente, creamos un objeto de tipo ts pasándole el argumento \( frequency=12 \) para indicar que se trata de una serie mensual:
precio <- read.csv("precioAzucarMex_2000_2013.csv")
precioMean <- cast(precio, anio + mes + tipo ~ X, value = "value", fun.aggregate = mean,
na.rm = TRUE)
ts <- ts(precioMean$DFCAIztapalapa[precioMean$tipo == "STANDAR"], start = c(2000,
10), frequency = 12)
Ahora, graficamos la serie:
plot(ts)
En la gráfica, claramente podemos descartar que se trate de una serie estacionaria, ya que el nivel no es constante. También se puede ver que aunque la serie tiene una tendencia creciente, ésta presenta varias intervenciones, por lo que se puede descartar tendencia lineal global. Para ver más claramente este comportamiento e investigar si la serie tiene componente estacional, graficamos la función de autocorrelación \( FAC \) y el diagrama de caja y brazos:
acf(ts)
boxplot(precioMean$DFCAIztapalapa ~ precioMean$mes)
En la \( FAC \) se puede notar que la serie tiene tendencia y no muestra estacionalidad; esto último se corrobora con el diagrama de caja y brazos.
Como se mencionó arriba, la serie ts muestra tendencia y no presenta estacionalidad. Esto sugiere que los posibles modelos a utilizar sean para series de nivel constante adaptivo, tendencia lineal estable o tendencia lineal adaptiva.
Para descartar alguno de los modelos anteriores, comparamos las varianzas de la serie original, una diferencia regular y 2 diferencias regulares:
var(ts)
## [1] 14957
dif_ts <- diff(ts, lag = 1, differences = 1)
var(dif_ts)
## [1] 783.3
var(diff(dif_ts))
## [1] 1326
La serie con la menor varianza es la de una diferencia regular, lo que sugiere que el modelo más adecuado es alguno para series de nivel constante adaptivo. Graficamos la serie de una diferencia regular y su \( FAC \) como apoyo al análisis de varianzas:
plot(dif_ts)
acf(dif_ts)
En estas gráficas se puede observar que una diferencia regular eliminó la tendencia de la serie; que las autocorrelaciones, excepto la de orden 5, son todas cero; y que la varianza no es constante. Lo anterior sugiere aplicar una transformación potencia para estabilizar la varianza; quizá de esta forma también se elimine la autocorrelación de orden 5.
Primero hacemos la prueba Goldfeld-Quandt para corroborar que la serie es heteroscedástica:
gqtest(dif_ts ~ 1)
##
## Goldfeld-Quandt test
##
## data: dif_ts ~ 1
## GQ = 6.885, df1 = 78, df2 = 77, p-value = 1.329e-15
Notamos que se rechaza la hipótesis de homoscedasticidad.
Ahora, buscamos la transformación potencia que mejor ayude a estabilizar la varianza, es decir, la que menor coeficiente de variación muestra entre las subseries:
subseries <- lapply(2001:2012, function(x) window(ts, start = x, end = c(x,
12)))
medias <- unlist(lapply(subseries, mean))
desviaciones <- unlist(lapply(subseries, sd))
lambda <- c(-3, -2.5, -2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2)
coefVar_lambda <- lapply(lambda, function(x) desviaciones/medias^(1 - x))
CV <- unlist(lapply(coefVar_lambda, sd))/unlist(lapply(coefVar_lambda, mean))
names(CV) <- lambda
CV
## -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5
## 0.8170 0.7140 0.6459 0.6227 0.6447 0.7005 0.7752 0.8566 0.9376 1.0144
## 2
## 1.0862
La tabla de arriba muestra el coeficiente de variación para distintos valores de \( \lambda \), escogemos \( \lambda=-1.5 \) para estabilizar nuestra serie y la graficamos, junto con su diferencia regular.
transf_ts <- -ts^(-1.5)
plot(transf_ts)
plot(diff(transf_ts))
Notamos al inicio de la serie de primeras diferencias una variación “grande”. Aplicamos nuevamente la prueba Goldfeld-Quandt para ver si la serie es homoscedástica en los extremos:
gqtest(transf_ts ~ 1)
##
## Goldfeld-Quandt test
##
## data: transf_ts ~ 1
## GQ = 1.146, df1 = 78, df2 = 78, p-value = 0.2743
Como \( p-value=0.2743 \), podemos aceptar la hipótesis nula (homoscedasticidad).
Escogemos el modelo de suavizamiento exponencial simple, que es equivalente a un Holt-Winters sin tendencia y sin estacionalidad (notar que esto no quiere decir que la \( \beta \) y \( \gamma \) son cero, si no que simplemente no entran en el modelo).
HW <- HoltWinters(transf_ts, beta = FALSE, gamma = FALSE)
HW
## Holt-Winters exponential smoothing without trend and without seasonal component.
##
## Call:
## HoltWinters(x = transf_ts, beta = FALSE, gamma = FALSE)
##
## Smoothing parameters:
## alpha: 1
## beta : FALSE
## gamma: FALSE
##
## Coefficients:
## [,1]
## a -0.0001416
plot(HW)
plot(residuals(HW))
El hecho de que \( \alpha=1 \) en este modelo quiere decir que para estimar el nivel, únicamente se está tomando en cuenta el valor anterior de la serie. En pocas palabras, la historia de la serie no nos está dando información sobre el comportamiento de la misma, únicamente el último valor observado. Por lo tanto, los residuales son las primeras diferencias. Se tiene un modelo del tipo: \[ W_t = W_{t-1} + \epsilon_t \] Ahora, si \( \epsilon_t \) son normales, con media cero, varianza constante y no autocorrelacionados; se podría obtener un intervalo de predicción para los pronósticos, en otro caso, los pronósticos serían puntuales.
Heteroscedasticidad. Hacemos la prueba \( \textit{Goldfeld-Quandt} \).
gqtest(residuals(HW) ~ 1)
##
## Goldfeld-Quandt test
##
## data: residuals(HW) ~ 1
## GQ = 0.5539, df1 = 78, df2 = 77, p-value = 0.995
Aunque la prueba no rechaza la hipótesis de no heteroscedasticidad, en la gráfica de residuales se puede notar que la varianza no es constante, de hecho, si cortamos la serie para eliminar el pico que se observa en 2001, la prueba tiene un resultado opuesto:
gqtest(window(residuals(HW), start = 2002) ~ 1)
##
## Goldfeld-Quandt test
##
## data: window(residuals(HW), start = 2002) ~ 1
## GQ = 1.939, df1 = 71, df2 = 70, p-value = 0.003052
Concluimos que los residuales de la serie no son homoscedásticos.
Autocorrelación en los residuales. Graficamos la \( FAC \) y hacemos la prueba Ljung-Box:
acf(residuals(HW))
lags <- seq(24)
LB.test <- as.data.frame(t(sapply(lags, function(lag) Box.test(residuals(HW),
lag = lag, type = "Ljung-Box"))))
LB.test <- LB.test[, c("method", "parameter", "statistic", "p.value")]
LB.test
## method parameter statistic p.value
## 1 Box-Ljung test 1 7.529 0.00607
## 2 Box-Ljung test 2 7.56 0.02282
## 3 Box-Ljung test 3 7.581 0.05552
## 4 Box-Ljung test 4 8.697 0.06914
## 5 Box-Ljung test 5 9.822 0.08044
## 6 Box-Ljung test 6 10.85 0.09301
## 7 Box-Ljung test 7 11.55 0.1163
## 8 Box-Ljung test 8 11.73 0.1639
## 9 Box-Ljung test 9 12.4 0.1918
## 10 Box-Ljung test 10 12.4 0.2593
## 11 Box-Ljung test 11 12.52 0.3259
## 12 Box-Ljung test 12 12.65 0.3949
## 13 Box-Ljung test 13 12.72 0.4693
## 14 Box-Ljung test 14 13 0.5266
## 15 Box-Ljung test 15 13.4 0.5717
## 16 Box-Ljung test 16 14.07 0.5934
## 17 Box-Ljung test 17 16.34 0.4996
## 18 Box-Ljung test 18 18.43 0.4277
## 19 Box-Ljung test 19 18.97 0.4586
## 20 Box-Ljung test 20 19.37 0.4977
## 21 Box-Ljung test 21 19.54 0.5504
## 22 Box-Ljung test 22 19.82 0.5946
## 23 Box-Ljung test 23 20.29 0.6242
## 24 Box-Ljung test 24 20.3 0.6798
Si vemos los valores p de la prueba, para autocorrelaciones \( k\geq10 \) Concluimos que los residuales no presentan autocorrelación.
Normalidad. Graficas el histograma de los residuales y usamos las pruebas Kolmogorov-Smirnoff y Shapiro-Wilk:
hist(residuals(HW), 30)
ks.test(residuals(HW)/sd(residuals(HW)), y = "pnorm")
##
## One-sample Kolmogorov-Smirnov test
##
## data: residuals(HW)/sd(residuals(HW))
## D = 0.1318, p-value = 0.008564
## alternative hypothesis: two-sided
shapiro.test(residuals(HW))
##
## Shapiro-Wilk normality test
##
## data: residuals(HW)
## W = 0.8957, p-value = 4.199e-09
Vemos que ambas pruebas rechazan la hipótesis de que los residuales se distribuyen normal.
Del análisis de residuales podemos concluir que no se cumplen los supuestos de ruido blanco, por lo que no será posible encontrar un intervalo de predicción.
Como se está ajustando un modelo de nivel constante adaptivo, donde todo el peso del suavizamiento recae en la observación anterior, tenemos que los pronósticos, están dados por. \[ \hat{W}_N(h)=W_N, \quad h={1,2,...} \] La gráfica con los pronósticos, en la variable transformada, con \( h=1,2,...,12 \) se muestra a continuación:
plot(HW, predicted.values = predict(HW, n.ahead = 10))
Ahora, obtenemos los residuales en la escala original y obtenemos el \( EM \), \( ECM \), \( EAM \) dinámicos con:
\( lag=1 \)
h <- 1
ventanaPrueba <- (length(ts)/2):(length(ts) - h)
residuales <- c()
for (i in ventanaPrueba) {
residuales <- c(residuales, ts[i] - ts[i + h])
}
sum(residuales)/length(residuales) # EM
## [1] -0.3463
sum(residuales^2)/length(residuales) # ECM
## [1] 1352
sum(abs(residuales))/length(residuales) # EAM
## [1] 25.35
\( lag=6 \)
h <- 6
ventanaPrueba <- (length(ts)/2):(length(ts) - h)
residuales <- c()
for (i in ventanaPrueba) {
residuales <- c(residuales, ts[i] - ts[i + h])
}
sum(residuales)/length(residuales) # EM
## [1] -2.425
sum(residuales^2)/length(residuales) # ECM
## [1] 11421
sum(abs(residuales))/length(residuales) # EAM
## [1] 78.5
Los errores anteriores nos dan una idea de cómo, entre más nos alejamos de la última observación, el pronóstico empeora.
En series de tiempo donde no se observa tendencia ni estacionalidad, los modelos que mejor ajustan la serie son aquellos con nivel constante adaptivo. Sin embargo, cuando el nivel sufre cambios repentinos a lo largo del tiempo, el peso de las observaciones alejadas en la estimación del nivel tendera a cero, mientras el de las observaciones más cercanas a uno. En el extremo, todo el peso para la estimación del nivel se concentra en la observación anterior, resultando en modelos como este. En otras palabras, la serie no contiene información, además de la última observación, para pronosticar valores futuros.
Este comportamiento es de esperarse en series de precios (mercados perfectos) ya que tienden a autoregularse; los cambios se deben a factores endógenos.
Con este modelo, el precio promedio del azúcar en la Central de Abastos del DF en diciembre será de: \( 368.0585+/-25.35 \)