El análisis que se presenta corresponde a la práctica del modulo 2 del diplomado de Modelos Económetricos Dinámicos en el ITAM, impartido por la maestra María Esperanza Sainz López. El análisis fue realizado por Emiliano Díaz y Humberto González.
Se uso una serie de tiempo de lluvia acumulada mensual en el distrito federal. La serie se construye a partir de las observaciones diarias de 10 estaciones meteorológicas en ubicadas en el distrito federal. Primero se tomo el promedio diario de las 10 estaciones y después se sumo la lluvia total en el mes.
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(reshape))
suppressPackageStartupMessages(library(nortest))
suppressPackageStartupMessages(library(moments))
suppressPackageStartupMessages(library(lmtest))
Cargamos la serie de tiempo. Aunque contamos con datos desde 1960 solo usaremos los datos de enero de 1990 a diciembre de 2010.
setwd("/Users/EDSP/Documents/Diplomado/modulo2")
lluviaMensualDF2 <- read.csv("serieLluviaDF.csv", colClasses = c("Date", "numeric",
"numeric", "numeric"))
lluviaMensualDF <- lluviaMensualDF2[which(lluviaMensualDF2$date >= as.Date("1990-01-01")),
]
Gráfiquemos la serie de tiempo para el periodo completo:
p <- ggplot(lluviaMensualDF)
p <- p + geom_line(aes(x = date, y = Lluvia))
p
Como es muy larga, conviene graficar un pedazo de la serie para mejor apreciar sus caracteristicas:
indx <- which(lluviaMensualDF$date >= as.Date("2005-01-01") & lluviaMensualDF$date <=
as.Date("2010-12-31"))
p <- ggplot(lluviaMensualDF[indx, ])
p <- p + geom_line(aes(x = date, y = Lluvia))
p
Obtenemos diagramas de caja y bigotes por año y mes para empezar a aislar el componente de tendencia y estacional respectivamente.
boxplot(lluviaMensualDF$Lluvia ~ lluviaMensualDF$Year)
boxplot(lluviaMensualDF$Lluvia ~ lluviaMensualDF$Month)
Como es de esperarse no se puede observar una tendencia en un periodo de 20 años. Además, se puede observar claramente lo que cualquier capitalino ya sabe: de junio a septiembre la cuidad se la pasa de tromba en tromba. Notar que no solo hay estacionalidad en el nivel de la serie sino en su varianza. En ausencia de modelos de volatilidad tendremos que encontrar una transformación que ayude a estabilizar la varianza. Eso lo intentaremos más adelante. Mientras tanto visualizamos la estacionalidad de los primeros dos momentos graficando la media y varianza por mes:
DE <- cast(lluviaMensualDF, Month ~ ., value = "Lluvia", fun.aggregate = sd)
colnames(DE) <- c("Month", "sigma")
Prom <- cast(lluviaMensualDF, Month ~ ., value = "Lluvia", fun.aggregate = mean)
colnames(Prom) <- c("Month", "mu")
df <- merge(Prom, DE, by = "Month")
df <- melt(df, id.vars = "Month")
p <- ggplot(df)
p <- p + geom_line(aes(x = Month, y = value, colour = variable))
p <- p + geom_point(aes(x = Month, y = value))
p
Graficamos la función de autocorrelación y realizamos la prueba Ljung-Box para verificar si hay tendencia y de que orden es la estacionalidad.
ACF <- acf(lluviaMensualDF$Lluvia)
plot(ACF)
lags <- seq(24)
LB.test <- as.data.frame(t(sapply(lags, function(lag) Box.test(lluviaMensualDF$Lluvia,
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 116.4 0
## 2 Box-Ljung test 2 140.4 0
## 3 Box-Ljung test 3 142.1 0
## 4 Box-Ljung test 4 184.3 0
## 5 Box-Ljung test 5 270.8 0
## 6 Box-Ljung test 6 371.8 0
## 7 Box-Ljung test 7 451.1 0
## 8 Box-Ljung test 8 488.5 0
## 9 Box-Ljung test 9 489.1 0
## 10 Box-Ljung test 10 517.8 0
## 11 Box-Ljung test 11 622.6 0
## 12 Box-Ljung test 12 774.3 0
## 13 Box-Ljung test 13 870.4 0
## 14 Box-Ljung test 14 891.7 0
## 15 Box-Ljung test 15 893 0
## 16 Box-Ljung test 16 928.8 0
## 17 Box-Ljung test 17 1009 0
## 18 Box-Ljung test 18 1102 0
## 19 Box-Ljung test 19 1176 0
## 20 Box-Ljung test 20 1209 0
## 21 Box-Ljung test 21 1210 0
## 22 Box-Ljung test 22 1229 0
## 23 Box-Ljung test 23 1324 0
## 24 Box-Ljung test 24 1465 0
Se rechaza la hipotesis nula que la serie no está correlacionada con sus primeros 24 desfases. Ahora procedemos a obtener diferencias de nuestra serie de tiempo para indagar sobre la técnica extrapolativa más adecuada para modelarla. Primero calculamos la desviación etándar de la serie para poder determinar si la diferencia aplicada simplifica la estructura de la serie:
sd(lluviaMensualDF$Lluvia)
## [1] 70.7
Probamos diferencias con desfase de 1,2,…,24 meses y orden de 1,2,..,5. Para cada combinación calculamos la desviación estándar de la serie. Se muestran las curvas de nivel de la desviación estándar para diferencias de distinto desfase y orden y se muestran las combinaciones que resultan en series con menor desviación estándar:
sensibilidad <- expand.grid(lag = seq(24), differences = seq(5))
sensibilidad$deDifs <- sapply(1:dim(sensibilidad)[1], function(i) sd(diff(lluviaMensualDF$Lluvia,
lag = sensibilidad$lag[i], differences = sensibilidad$differences[i])))
head(sensibilidad[order(sensibilidad$deDifs, decreasing = F), ])
## lag differences deDifs
## 12 12 1 45.05
## 24 24 1 45.60
## 1 1 1 56.77
## 11 11 1 57.99
## 23 23 1 59.20
## 13 13 1 60.37
p <- ggplot(sensibilidad)
p <- p + stat_summary2d(aes(x = lag, y = differences, z = deDifs), fun = function(x) mean(x),
bins = 5)
p <- p + geom_contour(aes(x = lag, y = differences, z = deDifs), colour = "white",
linemitr = 1)
p <- p + xlab("desfase de diferencias") + ylab("orden de diferncias")
p <- p + guides(fill = guide_legend(title = "desviación estándar"))
p
Pasamos de una desviación estándar para la serie original de 70.7035 a una desviación estándar para \( S(12,1) \) de 45.0466.
Probamos si podemos reducir la varianza aún más si aplicamos una segunda diferencia. Procedemos como antes.
diffs12 <- diff(lluviaMensualDF$Lluvia, lag = 12, differences = 1)
sensibilidad <- expand.grid(lag = seq(24), differences = seq(5))
sensibilidad$deDifs <- sapply(1:dim(sensibilidad)[1], function(i) sd(diff(diffs12,
lag = sensibilidad$lag[i], differences = sensibilidad$differences[i])))
head(sensibilidad[order(sensibilidad$deDifs, decreasing = F), ])
## lag differences deDifs
## 10 10 1 57.69
## 1 1 1 58.82
## 21 21 1 60.19
## 20 20 1 60.56
## 16 16 1 61.25
## 18 18 1 61.28
p <- ggplot(sensibilidad)
p <- p + stat_summary2d(aes(x = lag, y = differences, z = deDifs), fun = function(x) mean(x),
bins = 5)
p <- p + geom_contour(aes(x = lag, y = differences, z = deDifs), colour = "white",
linemitr = 1)
p <- p + xlab("desfase de diferencias") + ylab("orden de diferncias")
p <- p + guides(fill = guide_legend(title = "desviación estándar"))
p
Observar que todas las combinaciones probadas resultan en series con mayor desviación etándar por lo que concluimos que la serie únicamente tiene un componente estacional con periodicidad de 12 meses.
Dadas las caracteristicas de las series probaremos los siguientes modelos de pronóstico
El modelo Holt-Winters aditivo se puede representar de la siguiente forma:
\[ W_t = \mu_t + E_t + a_t, t=1,...,N \\ \mu_t = \mu_{t-1} + \beta_{t-1} + \alpha_1a_t\\ \beta_t = \beta_{t-1} + \alpha_1\alpha_2a_t\\ E_t = E_{t-S} + \alpha_3(1-\alpha_1)a_t \]
En donde el nivel, \( \mu \), la tendencia \( \beta \) y el factor estacional \( E \) de cada periodo se calculan como un promedio ponderado según las siguientes ecuaciones:
\[ \hat{\mu}_t = \alpha_1(W_t - \hat{E}_{t-S}) + (1-\alpha_1)(\hat{\mu}_{t-1}+\hat{\beta}_{t-1})\\ \hat{\beta}_t = \alpha_2(\mu_t - \mu_{t-1}) + (1-\alpha_2)\beta_{t-1}\\ \hat{E}_t = \alpha_3(W_t - \hat{\mu}_t)+(1-\alpha_3)E_{t-S} \]
Ajustamos el modelo y mostramos \( \alpha_1 \), \( \alpha_2 \) y \( \alpha_3 \):
x <- ts(lluviaMensualDF$Lluvia, start = c(1990, 1), frequency = 12)
HW <- HoltWinters(x, seasonal = "additive")
HW$alpha #alfa1
## alpha
## 0.02466
HW$beta #alfa2
## beta
## 0.1136
HW$gamma #alfa3
## gamma
## 0.2459
Observar que \( \alpha_1 \) y \( \alpha_2 \) son muy chicos por lo que el nivel y tendencia de la serie cambiaran de forma muy lenta. Grafiquemos los tres componentes de la serie (nivel, tendencia y estacionalidad) y el ajuste a la serie (ajuste es línea roja).Tomamos el subperiodo 2000-2010:
plot(window(HW$fitted, start = c(2000, 1), end = c(2010, 12)))
y <- cbind(x, HW$fitted[, 1])
y <- cbind(y, y[, 1] - y[, 2])
y <- na.omit(y)
plot(window(y[, 1:2], start = c(2000, 1), end = c(2010, 12)), plot.type = "single",
col = c("black", "red"))
Observando la escala de cada componente es claro que el componente estacional es el más importante. Ahora empezamos a analizar los residuales para determinar si cumplen los supuestos de normalidad (\( \epsilon_t \sim N(0,\sigma^2) \)), homocedasticidad y no auto-correlación. Si los cumplen podremos calcular intervalos de confianza para los pronósticos. Empezamos graficando la serie de tiempo e histograma de los residuales.
plot(window(y[, 3], start = c(1990, 1), end = c(2010, 12)))
## Warning: 'start' value not changed
hist(y[, 3], 30)
Ya es aparente que los residuales no tienen una distribución normal dado que tienen una kurtosis muy elevada:
summary(y[, 3]) #media, minimo, maximo, cuantiles
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -165.00 -15.50 -2.16 1.46 20.30 185.00
skewness(y[, 3]) #sesgo, para normal debe ser 0
## [1] 0.2087
kurtosis(y[, 3]) #kurtosis-pearson, para normal debe de ser 3
## [1] 6.638
geary(y[, 3]) #kurtosis-geary, para normal debe ser 0.80
## [1] 0.7144
# qqplot
qqnorm((y[, 3] - mean(y[, 3]))/sd(y[, 3]))
abline(a = 0, b = 1, col = "red")
Realizamos varias pruebas de normalidad a los residuales:
# pruebas de normalidad pruebas basadas en la función de distribución
ad.test(y[, 3]) #anderson darling
##
## Anderson-Darling normality test
##
## data: y[, 3]
## A = 2.448, p-value = 3.327e-06
cvm.test(y[, 3]) #cramer-von mises
##
## Cramer-von Mises normality test
##
## data: y[, 3]
## W = 0.4615, p-value = 6.266e-06
# pruebas basadas en la función de distribución empirica
ks.test(as.numeric(y[, 3]), y = "pnorm", mean = mean(y[, 3]), sd = sd(y[, 3])) #kolmogorov-smirnoff
##
## One-sample Kolmogorov-Smirnov test
##
## data: as.numeric(y[, 3])
## D = 0.0894, p-value = 0.04314
## alternative hypothesis: two-sided
lillie.test(y[, 3]) #Lilliefors
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: y[, 3]
## D = 0.0894, p-value = 8.081e-05
pearson.test(y[, 3]) #Pearson chi-square
##
## Pearson chi-square normality test
##
## data: y[, 3]
## P = 34.95, p-value = 0.0025
sf.test(y[, 3]) #Shapiro-Francia
##
## Shapiro-Francia normality test
##
## data: y[, 3]
## W = 0.9485, p-value = 8.889e-07
shapiro.test(y[, 3]) #Shapiro-Wilk
##
## Shapiro-Wilk normality test
##
## data: y[, 3]
## W = 0.9554, p-value = 9.159e-07
agostino.test(y[, 3]) #prueba D'agostino para sesgo
##
## D'Agostino skewness test
##
## data: y[, 3]
## skew = 0.2087, z = 0.8865, p-value = 0.3754
## alternative hypothesis: data have a skewness
anscombe.test(y[, 3]) #prueba Anscombe-Glynn para medida pearson de kurtosis normal (kurtosis=3)
##
## Anscombe-Glynn kurtosis test
##
## data: y[, 3]
## kurt = 6.638, z = 5.127, p-value = 2.937e-07
## alternative hypothesis: kurtosis is not equal to 3
bonett.test(y[, 3])
##
## Bonett-Seier test for Geary kurtosis
##
## data: y[, 3]
## tau = 26.524, z = 6.461, p-value = 1.041e-10
## alternative hypothesis: kurtosis is not equal to sqrt(2/pi)
jarque.test(as.numeric(y[, 3])) #prueba jarque bera para sesgo/simetria
##
## Jarque-Bera Normality Test
##
## data: as.numeric(y[, 3])
## JB = 134.1, p-value < 2.2e-16
## alternative hypothesis: greater
Casi todas las pruebas confirman que la distribución de los residuales no es normal por lo que solo podremos realizar pronosticos puntuales con este modelo. Realizamos pruebas de homocedasticidad (Goldfeld-Quandt) y auto-correlación (durbin-watson y Ljung-Box):
gqtest(y[, 3] ~ 1)
##
## Goldfeld-Quandt test
##
## data: y[, 3] ~ 1
## GQ = 0.5841, df1 = 119, df2 = 119, p-value = 0.9982
dwtest(y[, 3] ~ 1)
##
## Durbin-Watson test
##
## data: y[, 3] ~ 1
## DW = 1.599, p-value = 0.0009125
## alternative hypothesis: true autocorrelation is greater than 0
lags <- seq(24)
LB.test <- as.data.frame(t(sapply(lags, function(lag) Box.test(y[, 3], 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 9.753 0.00179
## 2 Box-Ljung test 2 9.961 0.00687
## 3 Box-Ljung test 3 9.995 0.01861
## 4 Box-Ljung test 4 10.59 0.03154
## 5 Box-Ljung test 5 11.82 0.03735
## 6 Box-Ljung test 6 13.46 0.03623
## 7 Box-Ljung test 7 13.72 0.05632
## 8 Box-Ljung test 8 14.88 0.06145
## 9 Box-Ljung test 9 15 0.09103
## 10 Box-Ljung test 10 17.87 0.05721
## 11 Box-Ljung test 11 18.33 0.07417
## 12 Box-Ljung test 12 18.43 0.1033
## 13 Box-Ljung test 13 18.44 0.1414
## 14 Box-Ljung test 14 19.32 0.153
## 15 Box-Ljung test 15 19.71 0.1835
## 16 Box-Ljung test 16 19.94 0.2228
## 17 Box-Ljung test 17 20.75 0.2377
## 18 Box-Ljung test 18 21.36 0.2618
## 19 Box-Ljung test 19 21.54 0.3076
## 20 Box-Ljung test 20 21.62 0.3614
## 21 Box-Ljung test 21 22.12 0.3928
## 22 Box-Ljung test 22 25.8 0.2602
## 23 Box-Ljung test 23 25.87 0.3071
## 24 Box-Ljung test 24 26.45 0.3305
ACF <- acf(y[, 3])
plot(ACF) #pareciera que todavía hay cierta correlación con el lag-1
No rechazamos hipotesis de residuales homocedasticos. Además aunque la prueba durbin watson indica que hay autocorrelación de orden 1 la prueba Ljung-Box no rechaza la hipotesis nula que la serie no está correlacionada con sus primeros 24 desfases.
Ahora procedemos con la simulación de pronósticos para el periodo 1995-2010 (usamos el periodo 1990-1995 para entrenar al modelo). Se calculan pronósticos a un mes, seis meses, 12 meses y 60 meses.
for (year in 1995:2010) {
for (month in 1:12) {
ventana <- window(x, start = c(1990, 1), end = c(year, month))
modelo <- HoltWinters(ventana)
pronosticos <- predict(modelo, n.ahead = 12 * 5)
# inicializamos pronosticos
if (month == 1 & year == 1995) {
pronosticos1 <- window(pronosticos, start = c(year, month + 1),
end = c(year, month + 1))
pronosticos6 <- window(pronosticos, start = c(year, month + 6),
end = c(year, month + 6))
pronosticos12 <- window(pronosticos, start = c(year, month + 12),
end = c(year, month + 12))
pronosticos60 <- window(pronosticos, start = c(year, month + 60),
end = c(year, month + 60))
} else {
pronosticos1 <- ts(c(pronosticos1, window(pronosticos, start = c(year,
month + 1), end = c(year, month + 1))), start = c(1995, 2),
frequency = 12)
pronosticos6 <- ts(c(pronosticos6, window(pronosticos, start = c(year,
month + 6), end = c(year, month + 6))), start = c(1995, 7),
frequency = 12)
pronosticos12 <- ts(c(pronosticos12, window(pronosticos, start = c(year,
month + 12), end = c(year, month + 12))), start = c(1995, 13),
frequency = 12)
pronosticos60 <- ts(c(pronosticos60, window(pronosticos, start = c(year,
month + 60), end = c(year, month + 60))), start = c(1995, 61),
frequency = 12)
}
}
}
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
y1 <- cbind(x, pronosticos1)
y1 <- cbind(y1, y1[, 1] - y1[, 2])
y1 <- na.omit(y1)
y6 <- cbind(x, pronosticos6)
y6 <- cbind(y6, y6[, 1] - y6[, 2])
y6 <- na.omit(y6)
y12 <- cbind(x, pronosticos12)
y12 <- cbind(y12, y12[, 1] - y12[, 2])
y12 <- na.omit(y12)
y60 <- cbind(x, pronosticos1)
y60 <- cbind(y60, y60[, 1] - y60[, 2])
y60 <- na.omit(y60)
Obtenemos el error cuadratico medio, error absoluto medio y gráfica de pronósticos (en rojo) para el prónostico adelantado de un mes:
mean(y1[, 3]^2) #error cuadratico medio 1 mes
## [1] 1348
mean(abs(y1[, 3])) #error absoluto medio 1 mes
## [1] 25.28
plot(window(y1[, 1:2], start = c(2000, 1), end = c(2010, 12)), plot.type = "single",
col = c("black", "red"), main = "pronóstico adelantado 1 mes")
Pronóstico adelantado de 6, 12 y 60 meses:
mean(y6[, 3]^2) #error cuadratico medio 6 meses
## [1] 1382
mean(abs(y6[, 3])) #error absoluto medio 6 meses
## [1] 26.01
plot(window(y6[, 1:2], start = c(2000, 1), end = c(2010, 12)), plot.type = "single",
col = c("black", "red"), main = "pronóstico adelantado 6 meses")
mean(y12[, 3]^2) #error cuadratico medio 12 meses
## [1] 1391
mean(abs(y12[, 3])) #error absoluto medio 12 meses
## [1] 26.16
plot(window(y12[, 1:2], start = c(2000, 1), end = c(2010, 12)), plot.type = "single",
col = c("black", "red"), main = "pronóstico adelantado 12 meses")
mean(y60[, 3]^2) #error cuadratico medio 60 meses
## [1] 1348
mean(abs(y60[, 3])) #error absoluto medio 60 meses
## [1] 25.28
plot(window(y60[, 1:2], start = c(2000, 1), end = c(2010, 12)), plot.type = "single",
col = c("black", "red"), main = "pronóstico adelantado 60 meses")
Notar que como el componente más importante es estacional los pronósticos muy adelantados son igualmente buenos que aquellos más inmediatos.
Dado que los errores no son un proceso de ruido blanco intentamos una transformación sobre la serie original para obtener residuales que cumplan supuestos de normalidad, homocedasticidad y no auto correlación. Primero intentamos transformaciones potencia de distinto tipo pero no se logro obtener residuales que cumplieran los supuestos. Los resultados se omiten para no extender demasiado la extensión del análisis. Como experimento probamos una transformación alternativa que dio buenos resultados.Construimos la función de distribución empírica de la serie de lluvia y se la aplicamos a la serie para obtener una serie con distribución uniforme. Primero definimos funciones que implementan la transformación (toUniform) y su inversa (toRain) que usaremos para convertir los pronósticos de vuelta a la escala original (en este caso no supimos como detectar y corregir sesgo):
toUniform <- function(vec, vecOriginal) {
p0 <- sum(vecOriginal == 0)/length(vecOriginal)
lower <- p0 #min(vec[which(vec>0)])
upper <- 1 - lower
f <- ecdf(vecOriginal[which(vecOriginal > 0)])
res <- f(vec)
summary(res)
res <- lower + upper * res
res[which(vec == 0)] <- runif(sum(vec == 0), min = 0, max = p0)
res[which(res == 1)] <- 0.99
return(res)
}
toRain <- function(vec, orig) {
p0 <- sum(orig == 0)/length(orig)
ps <- seq(0, 1, 0.005)
qs <- quantile(orig, probs = ps)
indx <- findInterval(vec, ps)
summary(indx)
indx[which(indx == 0)] <- 1
indx[which(indx > length(qs))] <- length(qs)
res <- (qs[indx] + qs[indx + 1])/2
res[which(vec < p0)] <- 0
res[which(vec >= 1)] <- max(orig)
return(res)
}
A continuación se presenta la gráfica de la serie de tiempo, el histograma mostrando la distribución empirica y un boxplot mostrando la distribución por mes. Para cada gráfica se muestra la serie original y luego la serie transformada para poder compararlas:
plot(lluviaMensualDF$Lluvia, type = "l", main = "serie original")
plot(toUniform(lluviaMensualDF$Lluvia, lluviaMensualDF2$Lluvia), type = "l",
main = "serie transformada")
hist(lluviaMensualDF$Lluvia, 100, main = "distribución original")
hist(toUniform(vec = lluviaMensualDF$Lluvia, lluviaMensualDF2$Lluvia), 100,
main = "distribución transformada")
boxplot(lluviaMensualDF$Lluvia ~ lluviaMensualDF$Month, main = "distribución original")
boxplot(toUniform(lluviaMensualDF$Lluvia, lluviaMensualDF2$Lluvia) ~ lluviaMensualDF$Month,
main = "distribución transformada")
Ajustamos el modelo Holt-Winters estacional aditivo a los datos transformados. Obtuvimos mejores resultados quitando el componente de tendencia del modelo (beta=FALSE).
x <- toUniform(lluviaMensualDF$Lluvia, lluviaMensualDF2$Lluvia)
x <- ts(x, start = c(1990, 1), frequency = 12)
HW <- HoltWinters(x, seasonal = "additive", beta = FALSE)
y <- cbind(x, HW$fitted[, 1])
y <- cbind(y, y[, 1] - y[, 2])
y <- na.omit(y)
Las \( \alpha_1 \), \( \alpha_2 \) y \( \alpha_3 \) ajustadas son las siguientes:
HW$alpha #alfa1
## alpha
## 0.005664
HW$beta #alfa2
## [1] FALSE
HW$gamma #alfa3
## gamma
## 0.2987
A continuación se muestra la descomposición de la serie y e ajuste del modelo.
plot(window(HW$fitted, start = c(1990, 1), end = c(2010, 12)))
## Warning: 'start' value not changed
plot(window(y[, 1:2], start = c(2000, 1), end = c(2010, 12)), plot.type = "single",
col = c("black", "red"))
Realizamos análisis de residuales:
plot(window(y[, 3], start = c(1991, 1), end = c(2010, 12)))
hist(y[, 3], 30)
# algunos estadísticos para indagar normalidad
summary(y[, 3])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.5130 -0.1000 -0.0020 -0.0051 0.0763 0.4870
skewness(y[, 3]) #sesgo, para normal debe ser 0
## [1] 0.02046
kurtosis(y[, 3]) #kurtosis-pearson, para normal debe de ser 3
## [1] 4.07
geary(y[, 3])
## [1] 0.7607
# qqplot
qqnorm((y[, 3] - mean(y[, 3]))/sd(y[, 3]))
abline(a = 0, b = 1, col = "red")
Realizamos algunas pruebas de normalidad:
# pruebas de normalidad
ad.test(y[, 3]) #anderson darling
##
## Anderson-Darling normality test
##
## data: y[, 3]
## A = 0.7114, p-value = 0.06261
cvm.test(y[, 3]) #cramer-von mises
##
## Cramer-von Mises normality test
##
## data: y[, 3]
## W = 0.1041, p-value = 0.09802
ks.test(as.numeric(y[, 3]), y = "pnorm", mean = mean(y[, 3]), sd = sd(y[, 3])) #ks
##
## One-sample Kolmogorov-Smirnov test
##
## data: as.numeric(y[, 3])
## D = 0.0572, p-value = 0.4114
## alternative hypothesis: two-sided
lillie.test(y[, 3]) #Lilliefors
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: y[, 3]
## D = 0.0572, p-value = 0.05469
pearson.test(y[, 3]) #Pearson chi-square
##
## Pearson chi-square normality test
##
## data: y[, 3]
## P = 19.2, p-value = 0.2048
sf.test(y[, 3]) #Shapiro-Francia
##
## Shapiro-Francia normality test
##
## data: y[, 3]
## W = 0.9859, p-value = 0.01929
shapiro.test(y[, 3]) #Shapiro-Wilk
##
## Shapiro-Wilk normality test
##
## data: y[, 3]
## W = 0.9882, p-value = 0.04643
agostino.test(y[, 3]) #prueba D'agostino para sesgo
##
## D'Agostino skewness test
##
## data: y[, 3]
## skew = 0.0205, z = 0.0878, p-value = 0.9301
## alternative hypothesis: data have a skewness
anscombe.test(y[, 3]) #prueba Anscombe-Glynn para medida pearson de kurtosis normal (kurtosis=3)
##
## Anscombe-Glynn kurtosis test
##
## data: y[, 3]
## kurt = 4.070, z = 2.631, p-value = 0.008517
## alternative hypothesis: kurtosis is not equal to 3
bonett.test(y[, 3])
##
## Bonett-Seier test for Geary kurtosis
##
## data: y[, 3]
## tau = 0.1104, z = 2.7943, p-value = 0.005201
## alternative hypothesis: kurtosis is not equal to sqrt(2/pi)
jarque.test(as.numeric(y[, 3])) #prueba jarque bera para sesgo/simetria
##
## Jarque-Bera Normality Test
##
## data: as.numeric(y[, 3])
## JB = 11.47, p-value = 0.003238
## alternative hypothesis: greater
Hay cierta evidencia para rechazar normalidad pero mucho menos que con residuos del modelo anterior por lo que intentaremos construir intervalos de probabilidad para los pronósticos y en la simulación de pronósticos podremos checar si se cumplen.
Se presentan pruebas para la homocedasticidad y auto-correlación de los resiudales.
gqtest(y[, 3] ~ 1)
##
## Goldfeld-Quandt test
##
## data: y[, 3] ~ 1
## GQ = 0.7297, df1 = 119, df2 = 119, p-value = 0.9565
dwtest(y[, 3] ~ 1)
##
## Durbin-Watson test
##
## data: y[, 3] ~ 1
## DW = 1.595, p-value = 0.0008092
## alternative hypothesis: true autocorrelation is greater than 0
lags <- seq(24)
LB.test <- as.data.frame(t(sapply(lags, function(lag) Box.test(y[, 3], 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 9.807 0.001738
## 2 Box-Ljung test 2 12.06 0.0024
## 3 Box-Ljung test 3 13.03 0.004564
## 4 Box-Ljung test 4 13.48 0.009158
## 5 Box-Ljung test 5 14.05 0.01531
## 6 Box-Ljung test 6 16.72 0.01039
## 7 Box-Ljung test 7 16.75 0.0191
## 8 Box-Ljung test 8 17.16 0.02852
## 9 Box-Ljung test 9 17.25 0.04487
## 10 Box-Ljung test 10 17.27 0.06862
## 11 Box-Ljung test 11 17.27 0.1001
## 12 Box-Ljung test 12 17.29 0.1391
## 13 Box-Ljung test 13 18.5 0.1394
## 14 Box-Ljung test 14 22.24 0.07388
## 15 Box-Ljung test 15 26.47 0.03341
## 16 Box-Ljung test 16 27.22 0.03906
## 17 Box-Ljung test 17 28.45 0.03993
## 18 Box-Ljung test 18 29.22 0.0458
## 19 Box-Ljung test 19 30.6 0.04465
## 20 Box-Ljung test 20 30.87 0.05696
## 21 Box-Ljung test 21 30.97 0.07415
## 22 Box-Ljung test 22 33.98 0.04935
## 23 Box-Ljung test 23 34.34 0.06042
## 24 Box-Ljung test 24 34.34 0.07876
ACF <- acf(y[, 3])
plot(ACF)
No hay evidencia de heterocedasticidad pero hay cierta evidencia de auto-correlación. Tendremos que tener cuidado para corroborar que intevalos de predicción sean precisos.
Realizamos la simulación de los pronósticos con adelanto de 1, 6, 12 y 60 meses como antes. Construimos intervalos de predicción al 98% de probabilidad. Una vez realizados los pronósticos en la escala transformada se convierten de regreso a la escala original.
for (year in 1995:2010) {
for (month in 1:12) {
ventana <- window(x, start = c(1990, 1), end = c(year, month))
modelo <- HoltWinters(ventana)
pronosticos <- predict(modelo, n.ahead = 12 * 5, prediction.interval = T,
level = 0.98)
# inicializamos pronosticos
if (month == 1 & year == 1995) {
pronosticos1 <- window(pronosticos, start = c(year, month + 1),
end = c(year, month + 1))
pronosticos6 <- window(pronosticos, start = c(year, month + 6),
end = c(year, month + 6))
pronosticos12 <- window(pronosticos, start = c(year, month + 12),
end = c(year, month + 12))
pronosticos60 <- window(pronosticos, start = c(year, month + 60),
end = c(year, month + 60))
} else {
pronosticos1 <- ts(rbind(pronosticos1, window(pronosticos, start = c(year,
month + 1), end = c(year, month + 1))), start = c(1995, 2),
frequency = 12)
pronosticos6 <- ts(rbind(pronosticos6, window(pronosticos, start = c(year,
month + 6), end = c(year, month + 6))), start = c(1995, 7),
frequency = 12)
pronosticos12 <- ts(rbind(pronosticos12, window(pronosticos, start = c(year,
month + 12), end = c(year, month + 12))), start = c(1995, 13),
frequency = 12)
pronosticos60 <- ts(rbind(pronosticos60, window(pronosticos, start = c(year,
month + 60), end = c(year, month + 60))), start = c(1995, 61),
frequency = 12)
}
}
}
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning: optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
# creamos series para comparar pronosticos con observaciones reales
y1.t <- cbind(x, pronosticos1)
y1.t <- cbind(y1.t, y1.t[, 1] - y1.t[, 2])
y1.t <- na.omit(y1.t)
z1 <- as.data.frame(y1.t)
z1[, 1] <- toRain(z1[, 1], lluviaMensualDF$Lluvia)
z1[, 2] <- toRain(z1[, 2], lluviaMensualDF$Lluvia)
z1[, 3] <- toRain(z1[, 3], lluviaMensualDF$Lluvia)
z1[, 4] <- toRain(vec = z1[, 4], orig = lluviaMensualDF$Lluvia)
z1[, 5] <- z1[, 1] - z1[, 2]
z1 <- ts(z1, start = start(y1.t), frequency = 12)
y6.t <- cbind(x, pronosticos6)
y6.t <- cbind(y6.t, y6.t[, 1] - y6.t[, 2])
y6.t <- na.omit(y6.t)
z6 <- as.data.frame(y6.t)
z6[, 1] <- toRain(z6[, 1], lluviaMensualDF$Lluvia)
z6[, 2] <- toRain(z6[, 2], lluviaMensualDF$Lluvia)
z6[, 3] <- toRain(z6[, 3], lluviaMensualDF$Lluvia)
z6[, 4] <- toRain(vec = z6[, 4], orig = lluviaMensualDF$Lluvia)
z6[, 5] <- z6[, 1] - z6[, 2]
z6 <- ts(z6, start = start(y6.t), frequency = 12)
y12.t <- cbind(x, pronosticos12)
y12.t <- cbind(y12.t, y12.t[, 1] - y12.t[, 2])
y12.t <- na.omit(y12.t)
z12 <- as.data.frame(y12.t)
z12[, 1] <- toRain(z12[, 1], lluviaMensualDF$Lluvia)
z12[, 2] <- toRain(z12[, 2], lluviaMensualDF$Lluvia)
z12[, 3] <- toRain(z12[, 3], lluviaMensualDF$Lluvia)
z12[, 4] <- toRain(vec = z12[, 4], orig = lluviaMensualDF$Lluvia)
z12[, 5] <- z12[, 1] - z12[, 2]
z12 <- ts(z12, start = start(y12.t), frequency = 12)
y60.t <- cbind(x, pronosticos60)
y60.t <- cbind(y60.t, y60.t[, 1] - y60.t[, 2])
y60.t <- na.omit(y60.t)
z60 <- as.data.frame(y60.t)
z60[, 1] <- toRain(z60[, 1], lluviaMensualDF$Lluvia)
z60[, 2] <- toRain(z60[, 2], lluviaMensualDF$Lluvia)
z60[, 3] <- toRain(z60[, 3], lluviaMensualDF$Lluvia)
z60[, 4] <- toRain(vec = z60[, 4], orig = lluviaMensualDF$Lluvia)
z60[, 5] <- z60[, 1] - z60[, 2]
z60 <- ts(z60, start = start(y60.t), frequency = 12)
A continuación se presentan el error cuadratico medio, el error absoluto medio y los pronósticos a un mes de adelanto.
mean(z1[, 5]^2) #error cuadratico medio 1 mes
## [1] 1453
mean(abs(z1[, 5])) #error absoluto medio 1 mes
## [1] 24.28
plot(window(z1[, 1:4], start = c(2000, 1), end = c(2010, 12)), plot.type = "single",
col = c("black", "red", "blue", "blue"))
Observar que el error absoluto medio es un poco menor al que se obtuvo sin transformar los datos (25.28). Verifiquemos si se cumple el intervalo de probabilidad teórico de 98%:
sum(z1[, 1] > z1[, 4] & z1[, 1] < z1[, 3])/dim(z1)[1] * 100
## [1] 94.76
La desviación de los residuales de la normaliad causa que no se cumpla fielmente el intervalo de probabilidad teórico aunque consideramos que se pueden usar los intervalos mientras se tome en cuenta que no cubren el 98% si no alrededor del 95% de probabilidad.
A continuación se presentan los pronósticos adelantados a 6, 12 y 60 meses.
mean(z6[, 5]^2) #error cuadratico medio 6 meses
## [1] 1449
mean(abs(z6[, 5])) #error absoluto medio 6 meses
## [1] 24.82
sum(z6[, 1] > z6[, 4] & z6[, 1] < z6[, 3])/dim(z6)[1] * 100
## [1] 93.55
plot(window(z6[, 1:4], start = c(2000, 1), end = c(2010, 12)), plot.type = "single",
col = c("black", "red", "blue", "blue"))
mean(z12[, 5]^2) #error cuadratico medio 12 meses
## [1] 1491
mean(abs(z12[, 5])) #error absoluto medio 12 meses
## [1] 25.29
sum(z12[, 1] > z12[, 4] & z12[, 1] < z12[, 3])/dim(z12)[1] * 100
## [1] 94.44
plot(window(z12[, 1:4], start = c(2000, 1), end = c(2010, 12)), plot.type = "single",
col = c("black", "red", "blue", "blue"))
mean(z60[, 5]^2) #error cuadratico medio 12 meses
## [1] 2773
mean(abs(z60[, 5])) #error absoluto medio 12 meses
## [1] 33.52
sum(z60[, 1] > z60[, 4] & z60[, 1] < z60[, 3])/dim(z60)[1] * 100
## [1] 96.21
plot(window(z60[, 1:4], start = c(2000, 1), end = c(2010, 12)), plot.type = "single",
col = c("black", "red", "blue", "blue"))
En general el tamaño real del intervalo de probabilidad esta entre 94 y 95%. Los pronósticos que se logran con este modelo son mejores a los del modelo anterior. Comparemos los errores absolutos medios de ambos modelos:
a <- c(mean(abs(y1[, 3])), mean(abs(y6[, 3])), mean(abs(y12[, 3])), mean(abs(y60[,
3])))
b <- c(mean(abs(z1[, 5])), mean(abs(z6[, 5])), mean(abs(z12[, 5])), mean(abs(z60[,
5])))
res <- data.frame(adelantoPronostico = c(1, 6, 12, 60), modeloAnterior = a,
modeloTransformacion = b)
res
## adelantoPronostico modeloAnterior modeloTransformacion
## 1 1 25.28 24.28
## 2 6 26.01 24.82
## 3 12 26.16 25.29
## 4 60 25.28 33.52
Concluimos el análisis mencionando que pronósticos de este tipo para la lluvia podrian ser útiles para alguna agencia gubernamental como aquella que se encarga del drenaje. Con este análisis podrian realizar una buena planeación de la capacidad que debe tener el sistema de drenaje para poder enfrentar escenarios extremos.