Series de tiempo
Luis Apáez
23/2/2022
Estimación de la tendencia: ajuste de la curva
Consideremos los siguientes datos, los cuales considera el número de nacidos vivos (en miles) en USA de 1948 a 1979
library(forecast)
library(nortest)
library(lmtest)# librería necesaria
library(astsa)
plot(birth)
puede notarse un cierto tipo de tendencia polnomial. Podemos notar la presencia de un ciclo que se repite cada 12 meses, donde podemos auxiliarnos de los correlogramas
tsdisplay(birth)
Asimismo, podemos ver los componentes de la serie de tiempo
dec <- decompose(birth, type="additive")
plot(dec)
la cual sigue un modelo aditivo. Continuando con lo visto la clase pasada, podemos ajustar un modelo de regresión que describa la serie, esto es, realizar un ajuste sobre la tendencia de la serie. Para ello, obtenemos los valores correspondientes para el eje horizontal (\(t\)) para tener la asociación \(t-f(t)\)
t = time(birth) - 1947 lo anterior nos servirá para ajustar la tendencia. Luego, lo siguiente nos ayudará a obtener el ciclo estacional:
M <- factor(cycle(birth))
M## [1] 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1
## [26] 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2
## [51] 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3
## [76] 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4
## [101] 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5
## [126] 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6
## [151] 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7
## [176] 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8
## [201] 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9
## [226] 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10
## [251] 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11
## [276] 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12
## [301] 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1
## [326] 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2
## [351] 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1
## Levels: 1 2 3 4 5 6 7 8 9 10 11 12la cual registra la correspondiente parte del periodo en vez de las observaciones.
Ahora bien, como se dijo al inicio, parace que la tendencia sigue un cierto comportamiento polinomial, por lo que
# ajustaremos a un polinomio de grado 5
t2 = t ** 2
t3 = t ** 3
t4 = t ** 4
t5 = t ** 5
reg <- lm( birth ~ 0 + t + t2 + t3 + t4 + t5 + M , na.action=NULL )
summary(reg)##
## Call:
## lm(formula = birth ~ 0 + t + t2 + t3 + t4 + t5 + M, na.action = NULL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.300 -6.507 -0.222 5.656 38.933
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## t -6.840e+00 2.976e+00 -2.298 0.0221 *
## t2 3.208e+00 5.221e-01 6.144 2.15e-09 ***
## t3 -2.765e-01 3.888e-02 -7.112 6.30e-12 ***
## t4 8.506e-03 1.284e-03 6.622 1.31e-10 ***
## t5 -8.774e-05 1.550e-05 -5.661 3.10e-08 ***
## M1 2.933e+02 5.676e+00 51.673 < 2e-16 ***
## M2 2.728e+02 5.672e+00 48.090 < 2e-16 ***
## M3 2.964e+02 5.683e+00 52.158 < 2e-16 ***
## M4 2.759e+02 5.693e+00 48.461 < 2e-16 ***
## M5 2.869e+02 5.703e+00 50.311 < 2e-16 ***
## M6 2.876e+02 5.713e+00 50.334 < 2e-16 ***
## M7 3.138e+02 5.723e+00 54.828 < 2e-16 ***
## M8 3.212e+02 5.732e+00 56.029 < 2e-16 ***
## M9 3.172e+02 5.741e+00 55.238 < 2e-16 ***
## M10 3.107e+02 5.751e+00 54.037 < 2e-16 ***
## M11 2.911e+02 5.759e+00 50.539 < 2e-16 ***
## M12 3.005e+02 5.768e+00 52.097 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.38 on 356 degrees of freedom
## Multiple R-squared: 0.999, Adjusted R-squared: 0.9989
## F-statistic: 1.993e+04 on 17 and 356 DF, p-value: < 2.2e-16donde recordamos que la variable categórica M nos auxilia para poder detectar el ciclo en la ST. Del summary() anterior podemos inducir que el mes com mayor número de nacimientos vivos es en Agosto (M8) pues es el que tiene el número mayor en la parte de la estimación de los coeficientes; el mes con menos nacimientos vivos es febrero (M2). Notemos además que el coeficiente de determinación es \(R^{2}=0.999\) por lo que el modelo propuesto es bastante bueno, es decir, la propuesta de que la tendencia sigue el comportamiento de un polinomio de grado 5 es bastante buena. Además notemos que
anova(reg)## Analysis of Variance Table
##
## Response: birth
## Df Sum Sq Mean Sq F value Pr(>F)
## t 1 26081421 26081421 242204.16 < 2.2e-16 ***
## t2 1 7437134 7437134 69064.67 < 2.2e-16 ***
## t3 1 1973425 1973425 18326.14 < 2.2e-16 ***
## t4 1 357574 357574 3320.59 < 2.2e-16 ***
## t5 1 233197 233197 2165.57 < 2.2e-16 ***
## M 12 396951 33079 307.19 < 2.2e-16 ***
## Residuals 356 38335 108
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1el error cuadrático medio es de 108. Graficamente podemos ver que
plot(birth, type = "o")
lines(fitted(reg), col = 2, lwd = 3)
el modelo ajustado es bastante próximo al modelo original de los datos. El resultado obtenido (la gráfica de color rojo) se obtiene pues en el ajuste del modelo agregamos los valores correspondientes a la t, con los que podemos ajustar la componente de la tendencia, y los valores correspondientes a la M, con los que podemos ajustar la componente del ciclo estacional. Ahora bien, si es sólo de nuestro interés estimar y después ver la componente de la tendencia, quitamos de la regresión a la variable M:
# ajuste sin la componente del ciclo
reg2 <- lm( birth ~ 0 + t + t2 + t3 + t4 + t5, na.action=NULL )
# Graficamos nuevamente
plot(birth, type = "o")
lines(fitted(reg2), col = 2, lwd = 3)
asimismo podemos visualizar sólo el ciclo sin tendencia:
# ajuste sin la componente de la tendencia
reg3 <- lm( birth ~ 0 + M, na.action=NULL )
# Graficamos nuevamente
plot(birth, type = "o")
lines(fitted(reg3), col = 2, lwd = 3)
Ahora, es preciso mencionar que ajustar la componente de la tendencia a un polinomio de grado 5 fue mera intuición, no obstante, tenemos herramientas que nos permiten saber si la estimación fue buena, o si un modelo es “mejor” que otro. Para ello nos auxiliamos del coeficiente de determinación y del error cuadrático medio. Por ejemplo, consideremos el siguiente ajuste
reg4 <- lm( birth ~ 0 + t + t2 + M, na.action=NULL )
summary(reg4)##
## Call:
## lm(formula = birth ~ 0 + t + t2 + M, na.action = NULL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -43.690 -13.672 0.053 12.401 47.269
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## t 6.71789 0.43504 15.44 <2e-16 ***
## t2 -0.26596 0.01281 -20.76 <2e-16 ***
## M1 290.45258 4.30469 67.47 <2e-16 ***
## M2 268.73565 4.36146 61.62 <2e-16 ***
## M3 292.47123 4.36517 67.00 <2e-16 ***
## M4 272.08147 4.36879 62.28 <2e-16 ***
## M5 283.24379 4.37233 64.78 <2e-16 ***
## M6 283.99044 4.37578 64.90 <2e-16 ***
## M7 310.32144 4.37914 70.86 <2e-16 ***
## M8 317.84968 4.38242 72.53 <2e-16 ***
## M9 313.96226 4.38561 71.59 <2e-16 ***
## M10 307.69143 4.38871 70.11 <2e-16 ***
## M11 288.16623 4.39173 65.62 <2e-16 ***
## M12 297.74151 4.39467 67.75 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.81 on 359 degrees of freedom
## Multiple R-squared: 0.9969, Adjusted R-squared: 0.9968
## F-statistic: 8196 on 14 and 359 DF, p-value: < 2.2e-16anova(reg4)## Analysis of Variance Table
##
## Response: birth
## Df Sum Sq Mean Sq F value Pr(>F)
## t 1 26081421 26081421 82209.13 < 2.2e-16 ***
## t2 1 7437134 7437134 23441.99 < 2.2e-16 ***
## M 12 2885587 240466 757.95 < 2.2e-16 ***
## Residuals 359 113895 317
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ajustando la componente de la tendencia a un polinomio de grado 2, tenemos que el coeficiente de determinación ha disminuido, asimismo el error cuadrático medio ha aumentado. Todo lo anterior respecto a la \(R^{2}\) y ECM del modelo con el ajuste al polinomio de grado 5. Como puede intuirse, al realizar el ajuste a un modelo respecto a un polinomio de gradp 3, la \(R^{2}\) aumentará y el ECM disminuirá. Veamos ahora lo que ocurre con un polinomio de grado 6:
t6 = t ** 6
reg5 <- lm( birth ~ 0 + t + t2 + t3 + t4 + t5 + t6 + M, na.action=NULL )
summary(reg5)##
## Call:
## lm(formula = birth ~ 0 + t + t2 + t3 + t4 + t5 + t6 + M, na.action = NULL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.950 -6.448 -0.500 5.536 36.506
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## t -2.414e+01 5.343e+00 -4.519 8.47e-06 ***
## t2 7.733e+00 1.277e+00 6.054 3.59e-09 ***
## t3 -7.899e-01 1.381e-01 -5.719 2.28e-08 ***
## t4 3.680e-02 7.425e-03 4.957 1.11e-06 ***
## t5 -8.339e-04 1.935e-04 -4.309 2.13e-05 ***
## t6 7.537e-06 1.949e-06 3.867 0.000131 ***
## M1 3.133e+02 7.595e+00 41.247 < 2e-16 ***
## M2 2.930e+02 7.634e+00 38.375 < 2e-16 ***
## M3 3.166e+02 7.646e+00 41.412 < 2e-16 ***
## M4 2.961e+02 7.656e+00 38.682 < 2e-16 ***
## M5 3.072e+02 7.665e+00 40.078 < 2e-16 ***
## M6 3.078e+02 7.673e+00 40.118 < 2e-16 ***
## M7 3.340e+02 7.681e+00 43.492 < 2e-16 ***
## M8 3.414e+02 7.687e+00 44.418 < 2e-16 ***
## M9 3.374e+02 7.692e+00 43.863 < 2e-16 ***
## M10 3.310e+02 7.697e+00 43.003 < 2e-16 ***
## M11 3.113e+02 7.701e+00 40.426 < 2e-16 ***
## M12 3.207e+02 7.704e+00 41.632 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.18 on 355 degrees of freedom
## Multiple R-squared: 0.999, Adjusted R-squared: 0.9989
## F-statistic: 1.956e+04 on 18 and 355 DF, p-value: < 2.2e-16anova(reg5)## Analysis of Variance Table
##
## Response: birth
## Df Sum Sq Mean Sq F value Pr(>F)
## t 1 26081421 26081421 251698.66 < 2.2e-16 ***
## t2 1 7437134 7437134 71772.03 < 2.2e-16 ***
## t3 1 1973425 1973425 19044.53 < 2.2e-16 ***
## t4 1 357574 357574 3450.76 < 2.2e-16 ***
## t5 1 233197 233197 2250.47 < 2.2e-16 ***
## t6 1 133005 133005 1283.57 < 2.2e-16 ***
## M 12 265495 22125 213.51 < 2.2e-16 ***
## Residuals 355 36786 104
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1vemos que la \(R^{2}\) se mantiene prácticamente igual y el ECM ha disminuido de 108 a 104. Gráficamente
# ajuste sin la componente del ciclo
reg6 <- lm( birth ~ 0 + t + t2 + t3 + t4 + t5 + t6, na.action=NULL )
# Graficamos nuevamente
plot(birth, type = "o")
lines(fitted(reg6), col = 2, lwd = 3)