Examen, Series de Tiempo
Eduardo Peña A01285443_A00833282 Ramón Islas_Guillermo Delgado A00835495_Piero Guillen A00835869
rm(list=ls())Examen final de Series de Tiempo
El examen final es una parte cuestionario (10 pts) y otra parte este rmd (10 pts)
Realiza el total de los pasos que se solicitan en cada punto o paso, los puntos que se hacen mención en cada uno es basado en 100, pero representan el 10 pts (50% del examen)
1) Primer paso (5 pts)
Lista de librerías requeridas: tidyverse, forecast, quantmod, AER, MASS, tseries, stats, car, lmtest, urca, aTSA, openxlsx, ggfortify
Verifica la existencia de las librerias a utilizar, para lo cual
debes utilizar la siguiente función para cada librería
if(!require(libreria)) install.packages(libreria)
Carga de librerías
2) Carga y revisión de información que se solicita (15 pts)
Carga de la información del csv con nombre Serie_Examen,
en una variable con nombre ‘mi_serie’, la cual se deberá cargar con los
encabezados y tomar en cuenta el tipo de separación que es ‘;’
mi_serie <- read.csv("Serie_Examen.csv", header = TRUE, sep = ";")Ahora deberás transformar los datos en una serie de tiempo mediante
la función ts, cargar el resultado en una variabkle con
nombre ‘mi_serie_co’, tomar solo la columna CO2 e indicar que debe
iniciar en Enero/1959, con una frecuencia mensual.
mi_serie_co<- ts(mi_serie$CO2, start = c(1959, 01) , frequency = 12)
mi_serie_co## Jan Feb Mar Apr May Jun Jul Aug Sep Oct
## 1959 333.13 332.09 331.10 329.14 327.36 327.29 328.23 329.55 330.62 331.40
## 1960 333.92 333.43 331.85 330.01 328.51 328.41 329.25 330.97 331.60 332.60
## 1961 334.68 334.17 332.96 330.80 328.98 328.57 330.20 331.58 332.67 333.17
## 1962 336.82 336.12 334.81 332.56 331.30 331.22 332.37 333.49 334.71 335.23
## 1963 337.95 338.00 336.37 334.47 332.46 332.29 333.76 334.80 336.00 336.63
## 1964 339.05 339.27 337.64 335.68 333.77 334.09 335.29 336.76 337.77 338.26
## 1965 341.47 341.31 339.41 337.74 336.07 336.07 337.22 338.38 339.32 340.41
## 1966 343.02 342.54 340.88 338.75 337.05 337.13 338.45 339.85 340.90 341.70
## 1967 344.28 343.42 342.02 339.97 337.84 338.00 339.20 340.63 341.41 342.68
## 1968 345.92 345.40 344.16 342.11 340.11 340.15 341.38 343.02 343.87 344.59
## 1969 347.38 346.78 344.96 342.71 340.86 341.13 342.84 344.32 344.88 345.62
## 1970 348.53 347.87 346.00 343.86 342.55 342.57 344.11 345.49 346.04 346.70
## 1971 349.93 349.26 347.44 345.55 344.21 343.67 345.09 346.27 347.33 347.82
## 1972 351.71 350.94 349.10 346.77 345.73
## Nov Dec
## 1959 331.87 333.18
## 1960 333.57 334.72
## 1961 334.86 336.07
## 1962 336.54 337.79
## 1963 337.93 338.95
## 1964 340.10 340.88
## 1965 341.69 342.51
## 1966 342.70 343.65
## 1967 343.04 345.27
## 1968 345.11 347.07
## 1969 347.23 347.62
## 1970 347.38 349.38
## 1971 349.29 350.91
## 1972Revisar el resultado, observando los primeros datos de la serie
mean(mi_serie_co)## [1] 339.1325var(mi_serie_co)## [1] 36.56678sd(mi_serie_co)## [1] 6.047047Revisar la información cargada mediante plot, describir la tendencia, si es o no estacionario.
plot(mi_serie_co, main = "CO2 a lo largo de los años", ylab = "CO2")
Comentarios La serie de tiempo presenta un comportamiento estacional mas no estacionario, con valores máximos en el segundo semestre del año y mínimos en el primero. Hay una tendencia creciente con el tiempo, con una vairanza constante, y no se una caminata aleatoria. A lo largo de los diferentes análisis vamos a ser capaces de poder explicar mas sobre esta serie de tiempo
3) Analizar la información (10 pts)
Utilice la función acf de Dicke-Fuller, explique lo que observa en la gráfica
mi_serie_co_diff <- diff(mi_serie_co)
mi_serie_co_adf_con_d <- adf.test(mi_serie_co_diff)## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 -5.16 0.01
## [2,] 1 -8.12 0.01
## [3,] 2 -9.60 0.01
## [4,] 3 -7.41 0.01
## [5,] 4 -7.74 0.01
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -5.16 0.01
## [2,] 1 -8.14 0.01
## [3,] 2 -9.68 0.01
## [4,] 3 -7.52 0.01
## [5,] 4 -7.89 0.01
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -5.14 0.01
## [2,] 1 -8.11 0.01
## [3,] 2 -9.64 0.01
## [4,] 3 -7.48 0.01
## [5,] 4 -7.86 0.01
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01mi_serie_co_adf_sin_d <- adf.test(mi_serie_co_diff)## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 -5.16 0.01
## [2,] 1 -8.12 0.01
## [3,] 2 -9.60 0.01
## [4,] 3 -7.41 0.01
## [5,] 4 -7.74 0.01
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -5.16 0.01
## [2,] 1 -8.14 0.01
## [3,] 2 -9.68 0.01
## [4,] 3 -7.52 0.01
## [5,] 4 -7.89 0.01
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -5.14 0.01
## [2,] 1 -8.11 0.01
## [3,] 2 -9.64 0.01
## [4,] 3 -7.48 0.01
## [5,] 4 -7.86 0.01
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01adf_test_result <- ur.df(mi_serie_co, lags = 12, type = "trend", selectlags = "AIC")
summary(adf_test_result)##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.10308 -0.16342 0.00389 0.19355 0.83443
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 77.47596 34.46440 2.248 0.02622 *
## z.lag.1 -0.23465 0.10490 -2.237 0.02696 *
## tt 0.02859 0.01278 2.237 0.02694 *
## z.diff.lag1 -0.00963 0.12307 -0.078 0.93775
## z.diff.lag2 0.06677 0.12198 0.547 0.58505
## z.diff.lag3 -0.19844 0.11097 -1.788 0.07600 .
## z.diff.lag4 -0.08844 0.10622 -0.833 0.40659
## z.diff.lag5 -0.10899 0.09970 -1.093 0.27630
## z.diff.lag6 -0.16439 0.09459 -1.738 0.08454 .
## z.diff.lag7 -0.14531 0.09032 -1.609 0.11001
## z.diff.lag8 -0.19450 0.08634 -2.253 0.02592 *
## z.diff.lag9 -0.22976 0.08184 -2.807 0.00575 **
## z.diff.lag10 -0.32166 0.07732 -4.160 5.67e-05 ***
## z.diff.lag11 0.06683 0.08129 0.822 0.41250
## z.diff.lag12 0.43342 0.07712 5.620 1.08e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3288 on 133 degrees of freedom
## Multiple R-squared: 0.9398, Adjusted R-squared: 0.9334
## F-statistic: 148.2 on 14 and 133 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -2.2368 9.2107 2.5036
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -3.99 -3.43 -3.13
## phi2 6.22 4.75 4.07
## phi3 8.43 6.49 5.47acf(mi_serie_co, main = "Funcion de autocorrelacion de la serie de tiempo CO2")
> Comentarios > Como se puede observar el valor de p en los dos
test (con diferenciación y sin diferenciación) dan un valor menor de
0.05, significa que hay evidencia suficiente para rechazar la hipótesis
nula de que la serie temporal tiene una raíz unitaria, entonces podemos
concluir que la serie es estacionaria
4) Analizás los componentes por separado (10 pts)
Utilizar la función stl, informe si respaldan o no las
conclusiones (comentarios) de los puntos anteriores
plot(stl(mi_serie_co, s.window = "periodic"))
# ?stl
mi_serie_co_ts <- ts(mi_serie_co, start=c(1959,01), frequency=12)
# PFE_ts
plot(stl(mi_serie_co_ts, s.window = "periodic"))
>Comentarios >La serie de tiempo tiene una varianza alta ya que
sus valores estan muy dispersos de la media, una estacionalidad anual,
con picos en los meses de verano y caídas en los meses de invierno, asi
como una tendencia creciente.
5) Modelos ARIMA (20 pts)
Definir los valores de p,d y q, mediante las funciones vistas en clase, excepto auto.arima
ndiffs(mi_serie_co) #determina el valor de d=1## [1] 1Y_d1 <- diff(mi_serie_co) #se carga la diferencia
adf.test(Y_d1) #se puede detrminar que si es estacionaria## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 -5.16 0.01
## [2,] 1 -8.12 0.01
## [3,] 2 -9.60 0.01
## [4,] 3 -7.41 0.01
## [5,] 4 -7.74 0.01
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -5.16 0.01
## [2,] 1 -8.14 0.01
## [3,] 2 -9.68 0.01
## [4,] 3 -7.52 0.01
## [5,] 4 -7.89 0.01
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -5.14 0.01
## [2,] 1 -8.11 0.01
## [3,] 2 -9.64 0.01
## [4,] 3 -7.48 0.01
## [5,] 4 -7.86 0.01
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01ndiffs(Y_d1) #se verifica si se requiere otra diferencia## [1] 0Y_d1_rl <- lm(Y_d1~time(Y_d1))
summary(Y_d1_rl) #Se puede observar que tiene p-value cercano a 1, la r2 ajustada es negativa##
## Call:
## lm(formula = Y_d1 ~ time(Y_d1))
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4198 -1.1677 0.4334 1.0528 2.1475
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.230910 51.243494 -0.063 0.950
## time(Y_d1) 0.001684 0.026069 0.065 0.949
##
## Residual standard error: 1.269 on 158 degrees of freedom
## Multiple R-squared: 2.64e-05, Adjusted R-squared: -0.006303
## F-statistic: 0.004171 on 1 and 158 DF, p-value: 0.9486plot(Y_d1, type = "l", main = "Primera diferencia de CO2", ylab = "Cantidad", xlab = "Año")
ndiffs(mi_serie_co) #determina el valor de d=1## [1] 1d <- ndiffs(mi_serie_co)
print(paste("El valor de 'd' es", d))## [1] "El valor de 'd' es 1"El valor de ‘d’ es
pacf_values <- Pacf(diff(mi_serie_co, differences = d), plot = FALSE, lag.max = 20)
d_aprox <- max(which(abs(pacf_values$acf)>1.96/sqrt(length(mi_serie_co))))
print(paste("El valor aproximado de d es:", d_aprox))## [1] "El valor aproximado de d es: 13"El valor apróximado de ‘p’ es
acf_values <- acf(diff(mi_serie_co, differences = 1), plot = FALSE, lag.max = 20)
q_aprox <- max(which(abs(pacf_values$acf)>1.96/sqrt(length(mi_serie_co))))
print(paste("El valor aproximado de q es:", q_aprox))## [1] "El valor aproximado de q es: 13"El valor apróximado de ‘q’ es
Si ‘d’>0 entonces haga acf y pacf según el valor de ‘d’
Mostrar acf, comente el significado de la gráficas
Y_diff <- diff(mi_serie_co)
acf(Y_diff)
plot(ts(Y_diff, frequency = 1))
– Comentario ACF: Con esta gráfica se puede entender mucho más sobre la
fuerte estacionalidad de la serie de tiempo, dando una fuerte
correlación entre los valores pasados. Debido a que los valores tienen
una caida rápida a 0, nos da a enteder a que no existe una dependencia
lineal. –
Mostrar pacf, comente el significado de la gráficas
pacf(Y_diff)
plot(ts(Y_diff, frequency = 1))
– Comentario El PACF de la serie de tiempo indica una fuerte correlación
positiva directa entre las observaciones consecutivas,
independientemente de los valores previos. Al igual que con el ACF, no
se ve una existente dependencia lineal. –
6) Verifica tu información (10 pts)
Verifica con un auto.arima, cuales son los valores de p,d y q, anotalos al final
auto.arima(mi_serie_co)## Series: mi_serie_co
## ARIMA(1,0,1)(2,1,2)[12] with drift
##
## Coefficients:
## ar1 ma1 sar1 sar2 sma1 sma2 drift
## 0.8880 -0.4252 -0.4862 -0.0029 -0.2211 -0.5112 0.1197
## s.e. 0.0529 0.0949 0.3324 0.1608 0.3334 0.3084 0.0029
##
## sigma^2 = 0.09842: log likelihood = -34.46
## AIC=84.92 AICc=85.95 BIC=108.96Comentarios El valor de ar1 indica que el valor actual de la serie de CO2 está autocorrelacionado positivamente con su valor anterior, el ma1 indica que existe una correlación negativa con la primera diferencia de la serie, los sar1 y 2 muestran la estacionalidad anual en la serie, el drift muestra que existe una tendencia creciente y el AIC muestra que es un valor relativamente bajo, lo cual signifca que es la versión con más ajuste.
7) Desarrolla del modelo (15 pts)
Desarrolla el modelo ARIMA con los datos del punto anterior, carga los resultados en una variable con nombre ‘mi_modelo’
mi_modelo <- arima(mi_serie_co, order=c(1,0,1), seasonal = list(order = c(2,1,2), period = 12))Revisa sus residuales, comparalos con el proceso estacionario (ruido blanco) que encontraste en el paso 2 y muestra su valor AIC.
resid2 <- mi_modelo$residuals^2
plot(mi_modelo$residuals)
plot(resid2)
acf(mi_modelo$residuals)
acf(resid2)
pacf(mi_modelo$residuals)
pacf(resid2)
Comentarios Según el modelo de autocorrelación, varios valores superan la zona de valores no signficativos de correlación, mientras que el resto de valores dentro de las lineas azules muestran una correlación significativa. La mayoria de los valores en los modelos de residuales muestran una correlación significativa ya que se encuentran dentro de las lineas azules.
8) Proyecciones (10 pts)
Realiza predicciones para los próximos 5 años (la frecuencia de la serie es mensual) y muestra estas gráficamente. Para esto puedes utilizar el comando ‘predict’ o ‘forecast’ estudiado en clase.
proyecciones = forecast::forecast(mi_modelo, h=60)
proyecciones## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jun 1972 345.4190 345.0328 345.8053 344.8283 346.0098
## Jul 1972 346.8200 346.3886 347.2514 346.1602 347.4798
## Aug 1972 348.1627 347.6904 348.6349 347.4405 348.8849
## Sep 1972 349.0322 348.5224 349.5419 348.2525 349.8118
## Oct 1972 349.7497 349.2049 350.2944 348.9166 350.5828
## Nov 1972 350.9511 350.3735 351.5286 350.0678 351.8344
## Dec 1972 352.2753 351.6667 352.8839 351.3446 353.2061
## Jan 1973 352.9395 352.3015 353.5776 351.9637 353.9154
## Feb 1973 352.3244 351.6581 352.9906 351.3054 353.3433
## Mar 1973 350.6199 349.9266 351.3131 349.5596 351.6801
## Apr 1973 348.4795 347.7602 349.1987 347.3795 349.5795
## May 1973 347.0048 346.2605 347.7491 345.8665 348.1431
## Jun 1973 346.9089 346.1051 347.7128 345.6796 348.1383
## Jul 1973 348.2975 347.4565 349.1384 347.0114 349.5836
## Aug 1973 349.6737 348.7972 350.5502 348.3332 351.0141
## Sep 1973 350.4974 349.5867 351.4080 349.1047 351.8900
## Oct 1973 351.2209 350.2774 352.1644 349.7779 352.6638
## Nov 1973 352.2630 351.2878 353.2383 350.7715 353.7545
## Dec 1973 353.7178 352.7118 354.7238 352.1793 355.2563
## Jan 1974 354.3092 353.2734 355.3450 352.7251 355.8933
## Feb 1974 353.7069 352.6422 354.7716 352.0785 355.3352
## Mar 1974 352.0046 350.9117 353.0975 350.3332 353.6760
## Apr 1974 349.9240 348.8037 351.0444 348.2107 351.6374
## May 1974 348.3767 347.2297 349.5238 346.6225 350.1310
## Jun 1974 348.2547 347.0672 349.4421 346.4387 350.0706
## Jul 1974 349.6500 348.4308 350.8691 347.7854 351.5145
## Aug 1974 350.9990 349.7489 352.2490 349.0872 352.9108
## Sep 1974 351.8563 350.5761 353.1365 349.8985 353.8142
## Oct 1974 352.5631 351.2535 353.8727 350.5602 354.5660
## Nov 1974 353.7017 352.3633 355.0401 351.6548 355.7486
## Dec 1974 355.1056 353.7390 356.4721 353.0156 357.1955
## Jan 1975 355.7414 354.3473 357.1355 353.6093 357.8735
## Feb 1975 355.1230 353.7019 356.5441 352.9496 357.2964
## Mar 1975 353.4111 351.9636 354.8587 351.1973 355.6250
## Apr 1975 351.2884 349.8148 352.7619 349.0348 353.5420
## May 1975 349.8024 348.3034 351.3015 347.5098 352.0951
## Jun 1975 349.6809 348.1349 351.2269 347.3165 352.0453
## Jul 1975 351.0727 349.4927 352.6528 348.6562 353.4892
## Aug 1975 352.4331 350.8197 354.0465 349.9657 354.9006
## Sep 1975 353.2751 351.6291 354.9211 350.7577 355.7925
## Oct 1975 353.9894 352.3114 355.6674 351.4232 356.5557
## Nov 1975 355.0868 353.3775 356.7962 352.4726 357.7011
## Dec 1975 356.5087 354.7686 358.2489 353.8474 359.1700
## Jan 1976 357.1252 355.3549 358.8956 354.4177 359.8328
## Feb 1976 356.5136 354.7136 358.3137 353.7607 359.2666
## Mar 1976 354.8059 352.9767 356.6351 352.0083 357.6035
## Apr 1976 352.7007 350.8428 354.5587 349.8593 355.5422
## May 1976 351.1867 349.3005 353.0729 348.3020 354.0714
## Jun 1976 351.0657 349.1356 352.9957 348.1139 354.0174
## Jul 1976 352.4582 350.4936 354.4228 349.4536 355.4628
## Aug 1976 353.8136 351.8151 355.8122 350.7571 356.8701
## Sep 1976 354.6609 352.6290 356.6928 351.5534 357.7684
## Oct 1976 355.3716 353.3070 357.4363 352.2140 358.5293
## Nov 1976 356.4840 354.3871 358.5809 353.2770 359.6910
## Dec 1976 357.8989 355.7702 360.0275 354.6433 361.1544
## Jan 1977 358.5221 356.3621 360.6820 355.2187 361.8254
## Feb 1977 357.9072 355.7165 360.0979 354.5568 361.2576
## Mar 1977 356.1972 353.9761 358.4183 352.8004 359.5940
## Apr 1977 354.0847 351.8337 356.3357 350.6421 357.5273
## May 1977 352.5809 350.3005 354.8614 349.0933 356.0686plot(proyecciones)