Cargamos los datos:
suppressPackageStartupMessages(library(gdata))
suppressPackageStartupMessages(library(forecast))
data <- read.xls("Autos.xls")
data <- data[1:309, ]
serie <- ts(data[, 2], start = c(1988, 1), frequency = 12)
Graficamos la serie:
plot(serie)
La serie presenta tendencia global, sin embargo se observa un comportamiento local que parecería adaptivo. También se nota estacionalidad, sobre todo a partir del 2000 y que la varianza aumenta con el nivel.
Cortamos la gráfica para incluir únicamente datos de Abril de 1995 a junio de 2013.
serieW <- window(serie, start = c(1995, 4), end = c(2013, 6))
Ajustamos un modelo de regresión lineal para eliminar la tendencia de la serie. El modelo ajustado es \( y=\beta_0+\beta_1t+\beta_2t^2+\beta_3t^3+\beta_4t^4 \)
serieDF <- data.frame(t = 1:length(serieW), obs = as.numeric(serieW))
serieDF$t2 <- serieDF$t^2
serieDF$t3 <- serieDF$t^3
serieDF$t4 <- serieDF$t^4
serieW_lm <- lm(serieDF$obs ~ serieDF$t + serieDF$t2 + serieDF$t3 + serieDF$t4)
summary(serieW_lm)
##
## Call:
## lm(formula = serieDF$obs ~ serieDF$t + serieDF$t2 + serieDF$t3 +
## serieDF$t4)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15402 -4332 -1351 2169 39348
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.39e+03 2.78e+03 2.65 0.00856 **
## serieDF$t 4.10e+02 1.75e+02 2.35 0.01984 *
## serieDF$t2 1.08e+01 3.22e+00 3.35 0.00094 ***
## serieDF$t3 -1.32e-01 2.19e-02 -6.00 8.2e-09 ***
## serieDF$t4 3.62e-04 4.95e-05 7.33 4.8e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8020 on 214 degrees of freedom
## Multiple R-squared: 0.765, Adjusted R-squared: 0.76
## F-statistic: 174 on 4 and 214 DF, p-value: <2e-16
Vemos que se rechaza la hipótesis de que \( \beta_i=0 \forall i \) y las pruebas individuales sobre cada \( \beta_i \).
Graficamos la serie original contra \( \^y \) estimada en el punto anterior.
plot(serieDF$t, serieDF$obs, type = "l")
lines(serieDF$t, serieW_lm$fitted, col = "blue")
Graficamos también la serie sin tendencia, que no son más que los residuos del modelo anterior:
plot(serieW_lm$residuals, type = "l")
Ahora, graficamos la función de autocorrelación \( FAC \) de los residuos y obtenemos:
acf(serieW_lm$residuals)
Notamos estacionalidad de orden 12 en la serie sin tendencia.
Ahora, modelamos únicamente la estacionalidad sin eliminar la tendencia de la serie. Ajustamos el modelo \( y_{t}=\sum_{k=0}^{11}\alpha_{tk}x_{tk} \), donde \( x_{tk}=1 \) cuando \( modulo(t)=k \), es decir, cuando la onbservación \( t \) pertenece al mes \( k \) del año.
serieDF$mes <- (serieDF$t + 3)%%12
serieDF$mes[serieDF$mes == 0] <- 12
serieW_lm_est <- lm(serieDF$obs ~ as.factor(serieDF$mes) - 1)
summary(serieW_lm_est)
##
## Call:
## lm(formula = serieDF$obs ~ as.factor(serieDF$mes) - 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -52459 -8714 6079 10702 32891
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## as.factor(serieDF$mes)1 47918 3610 13.3 <2e-16 ***
## as.factor(serieDF$mes)2 44620 3610 12.4 <2e-16 ***
## as.factor(serieDF$mes)3 46233 3610 12.8 <2e-16 ***
## as.factor(serieDF$mes)4 38904 3514 11.1 <2e-16 ***
## as.factor(serieDF$mes)5 41440 3514 11.8 <2e-16 ***
## as.factor(serieDF$mes)6 42027 3514 12.0 <2e-16 ***
## as.factor(serieDF$mes)7 42554 3610 11.8 <2e-16 ***
## as.factor(serieDF$mes)8 44346 3610 12.3 <2e-16 ***
## as.factor(serieDF$mes)9 42059 3610 11.7 <2e-16 ***
## as.factor(serieDF$mes)10 45173 3610 12.5 <2e-16 ***
## as.factor(serieDF$mes)11 48004 3610 13.3 <2e-16 ***
## as.factor(serieDF$mes)12 66416 3610 18.4 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15300 on 207 degrees of freedom
## Multiple R-squared: 0.906, Adjusted R-squared: 0.901
## F-statistic: 166 on 12 and 207 DF, p-value: <2e-16
Graficamos el modelo ajustado (rojo) contra la serie original.
plot(serieDF$t, serieDF$obs, type = "l")
lines(serieDF$t, serieW_lm_est$fitted, col = "red")
Graficamos los residuos y la \( FAC \):
plot(serieW_lm_est$residuals, type = "l")
acf(serieW_lm_est$residuals)
Notamos que la tendencia se puede observar más fácilmente en los residuos que en la serie original. Así mismo, la gráfica de la \( FAC \) muestra la presencia de tendencia en los residuos.
Proponemos modelar la estacionalidad incluyendo variables indicadoras para cada mes en el modelo anterior:
serieW_lm_tend_est <- lm(serieDF$obs ~ serieDF$t + serieDF$t2 + serieDF$t3 +
serieDF$t4 + as.factor(serieDF$mes))
summary(serieW_lm_tend_est)
##
## Call:
## lm(formula = serieDF$obs ~ serieDF$t + serieDF$t2 + serieDF$t3 +
## serieDF$t4 + as.factor(serieDF$mes))
##
## Residuals:
## Min 1Q Median 3Q Max
## -18910 -2775 484 2816 18929
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.03e+04 1.98e+03 5.19 5.0e-07 ***
## serieDF$t 3.27e+02 1.03e+02 3.17 0.00175 **
## serieDF$t2 1.22e+01 1.90e+00 6.44 8.4e-10 ***
## serieDF$t3 -1.41e-01 1.29e-02 -10.90 < 2e-16 ***
## serieDF$t4 3.83e-04 2.92e-05 13.12 < 2e-16 ***
## as.factor(serieDF$mes)2 -3.56e+03 1.57e+03 -2.27 0.02449 *
## as.factor(serieDF$mes)3 -2.22e+03 1.57e+03 -1.41 0.15938
## as.factor(serieDF$mes)4 -7.82e+03 1.55e+03 -5.03 1.1e-06 ***
## as.factor(serieDF$mes)5 -5.57e+03 1.55e+03 -3.59 0.00042 ***
## as.factor(serieDF$mes)6 -5.27e+03 1.55e+03 -3.40 0.00082 ***
## as.factor(serieDF$mes)7 -3.88e+03 1.57e+03 -2.46 0.01454 *
## as.factor(serieDF$mes)8 -2.32e+03 1.57e+03 -1.48 0.14158
## as.factor(serieDF$mes)9 -4.85e+03 1.57e+03 -3.08 0.00233 **
## as.factor(serieDF$mes)10 -1.98e+03 1.57e+03 -1.26 0.20917
## as.factor(serieDF$mes)11 6.00e+02 1.57e+03 0.38 0.70315
## as.factor(serieDF$mes)12 1.88e+04 1.57e+03 11.93 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4720 on 203 degrees of freedom
## Multiple R-squared: 0.923, Adjusted R-squared: 0.917
## F-statistic: 162 on 15 and 203 DF, p-value: <2e-16
Graficamos la serie original contra el modelo ajustado que incluye tendencia y estacionalidad.
plot(serieDF$t, serieDF$obs, type = "l")
lines(serieDF$t, serieW_lm_tend_est$fitted, col = "red")
Ahora graficamos los residuos, los residuos al cuadrado y su \( FAC \):
plot(serieW_lm_tend_est$residuals, type = "l")
plot(serieW_lm_tend_est$residuals^2, type = "l")
acf(serieW_lm_tend_est$residuals)
En las gráficas de los residuos se puede notar que, en gran medida, el modelo ajusta tendencia y estacionalidad. Sin embargo, los residuos presentan heteroscedasticidad, loq ue quiere decir que la variabilidad no está siendo explicada completamente por el modelo. Más aún, es posible notar ciertas intervenciones en la gráfica de residuos al cuadrado. En conclusión, los residuos no son ruido blanco.
Ahora, hacemos el pronóstico para los próximos 3 meses:
predict1 <- serieW_lm_tend_est$coefficients[1] + serieW_lm_tend_est$coefficients[2] *
220 + serieW_lm_tend_est$coefficients[3] * 220^2 + serieW_lm_tend_est$coefficients[4] *
220^3 + serieW_lm_tend_est$coefficients[5] * 220^4 - 3877.155
predict2 <- serieW_lm_tend_est$coefficients[1] + serieW_lm_tend_est$coefficients[2] *
221 + serieW_lm_tend_est$coefficients[3] * 221^2 + serieW_lm_tend_est$coefficients[4] *
221^3 + serieW_lm_tend_est$coefficients[5] * 221^4 - 2320.866
predict3 <- serieW_lm_tend_est$coefficients[1] + serieW_lm_tend_est$coefficients[2] *
222 + serieW_lm_tend_est$coefficients[3] * 222^2 + serieW_lm_tend_est$coefficients[4] *
222^3 + serieW_lm_tend_est$coefficients[5] * 222^4 - 4849.714
Los pronósticos son \( jul2013 = 64,780.52 \), \( ago2013 = 67,899.22 \) y \( sep2013 = 66,994.62 \) y los valores observados para estos meses fueron de \( 56,640 \), \( 57,079 \), \( 50,732 \), respectivamente. La gráfica con los valores ajustados (rojo), pronosticados (azul) y los observados (negro) muestra a continuación:
plot(serieDF$t, serieDF$obs, type = "l")
lines(serieDF$t, serieW_lm_tend_est$fitted, col = "red")
lines(c(220, 221, 222), c(predict1, predict2, predict3), col = "blue")
lines(c(220, 221, 222), c(56640, 57079, 50732))
abline(v = 219)
De donde se ve que el error de pronóstico es alto: la tendencia en los valores observados es opuesta a los pronosticados.