rm(list=ls())
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)
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)
if(!require(tidyverse)) install.packages(tidyverse)
## Loading required package: tidyverse
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.4.4 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
if(!require(forecast)) install.packages(forecast)
## Loading required package: forecast
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
if(!require(quantmod)) install.packages(quantmod)
## Loading required package: quantmod
## Loading required package: xts
## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
##
## ######################### Warning from 'xts' package ##########################
## # #
## # The dplyr lag() function breaks how base R's lag() function is supposed to #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
## # source() into this session won't work correctly. #
## # #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
## # dplyr from breaking base R's lag() function. #
## # #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
## # #
## ###############################################################################
##
## Attaching package: 'xts'
##
## The following objects are masked from 'package:dplyr':
##
## first, last
##
## Loading required package: TTR
if(!require(AER)) install.packages(AER)
## Loading required package: AER
## Loading required package: car
## Loading required package: carData
##
## Attaching package: 'car'
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## some
##
## Loading required package: lmtest
## Loading required package: sandwich
## Loading required package: survival
if(!require(MASS)) install.packages(MASS)
## Loading required package: MASS
##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:dplyr':
##
## select
if(!require(tseries)) install.packages(tseries)
## Loading required package: tseries
if(!require(stats)) install.packages(stats)
if(!require(car)) install.packages(car)
if(!require(lmtest)) install.packages(lmtest)
if(!require(urca)) install.packages(urca)
## Loading required package: urca
if(!require(aTSA)) install.packages(aTSA)
## Loading required package: aTSA
##
## Attaching package: 'aTSA'
##
## The following objects are masked from 'package:tseries':
##
## adf.test, kpss.test, pp.test
##
## The following object is masked from 'package:forecast':
##
## forecast
##
## The following object is masked from 'package:graphics':
##
## identify
if(!require(openxlsx)) install.packages(openxlsx)
## Loading required package: openxlsx
if(!require(ggfortify)) install.packages(ggfortify)
## Loading required package: ggfortify
## Registered S3 methods overwritten by 'ggfortify':
## method from
## autoplot.Arima forecast
## autoplot.acf forecast
## autoplot.ar forecast
## autoplot.bats forecast
## autoplot.decomposed.ts forecast
## autoplot.ets forecast
## autoplot.forecast forecast
## autoplot.stl forecast
## autoplot.ts forecast
## fitted.ar forecast
## fortify.ts forecast
## residuals.ar forecast
Carga de librerÃas
tidyverse, forecast, quantmod, AER, MASS, tseries, stats, car, lmtest, urca, aTSA, openxlsx, ggfortify
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 ‘;’
Serie_Examen <- 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 variable con
nombre ‘mi_serie_co’, tomar solo la columna CO2 e indicar que debe
iniciar en Enero/1959, con una frecuencia mensual.
mi_serie <- read.csv('Serie_Examen.csv', header = TRUE, sep = ";")
mi_serie_co <- ts(mi_serie$CO2, start = c(1959, 01), frequency = 12)
print(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
## 1972
Revisar el resultado, observando los primeros datos de la serie
mean(mi_serie_co)
## [1] 339.1325
var(mi_serie_co)
## [1] 36.56678
sd(mi_serie_co)
## [1] 6.047047
Revisar la información cargada mediante plot, describir la tendencia, si es o no estacionario.
plot(mi_serie_co, main = "Serie de tiempo CO2", xlab = "Fecha", ylab = "CO2")
adf_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.47
#En este caso, si es estacionario ya que está incrementando y decreciendo gradualmente conforme pasan los años. Se puede observar una tendencia a la alza en esta serie de tiempo comenzando desde el año 1959 y terminando en 1972
Comentarios
Utilice la función acf de Dicke-Fuller, explique lo que observa en la gráfica
adf_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.47
acf(mi_serie_co, main = "Función de autocorrelación de la serie de tiempo CO2")
#Esta gráfica es una prueba estadÃstica que se usa para determinar si la serie tiene raÃz unitaria. Con esto podemos evaluar si esta serie tiene un movimiento estacionario o no. Después de analizar la gráfica, se muestra la estadÃstica de prueba en el eje vertical y el tamaño de la muestra en el eje horizontal. Esto muestra que tiene una tendencia baja, y que no es aleatorio ya que va por segmentos y se mueve con coordinación.
Comentarios
Utilizar la función stl, informe si respaldan o no las
conclusiones (comentarios) de los puntos anteriores
mi_serie <- read.csv("Serie_Examen.csv", header = TRUE, sep = ";")
mi_serie_co <- ts(mi_serie$CO2, start = c(1959, 1), frequency = 12)
head(mi_serie_co)
## Jan Feb Mar Apr May Jun
## 1959 333.13 332.09 331.10 329.14 327.36 327.29
plot(mi_serie_co, main = "Serie de tiempo CO2", xlab = "Año", ylab = "CO2")
adf.test(mi_serie_co)
## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 0.764 0.862
## [2,] 1 0.290 0.726
## [3,] 2 0.674 0.836
## [4,] 3 1.160 0.934
## [5,] 4 1.198 0.939
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -1.28 0.597
## [2,] 1 -3.50 0.010
## [3,] 2 -1.98 0.336
## [4,] 3 -1.51 0.516
## [5,] 4 -1.67 0.457
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -3.88 0.0168
## [2,] 1 -13.11 0.0100
## [3,] 2 -9.82 0.0100
## [4,] 3 -8.15 0.0100
## [5,] 4 -11.10 0.0100
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01
acf(mi_serie_co, main = "Función de autocorrelación")
#En esta parte el código analiza una serie de tiempo CO2, en donde la serie temporal muestra una tendencia a la alza. La otra gráfica, siendo la de autocorrelación muestra una disminución gradual siendo una correlación negativa.
Definir los valores de p,d y q, mediante las funciones vistas en clase, excepto auto.arima
ndiffs(mi_serie_co)
## [1] 1
d <- 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)
p_aprox <- max(which(abs(pacf_values$acf) > 1.96/sqrt(length(mi_serie_co))))
print(paste("El valor apróximado de 'p' es", p_aprox))
## [1] "El valor apróximado de 'p' 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(acf_values$acf) > 1.96/sqrt(length(mi_serie_co))))
print(paste("El valor apróximado de 'q' es", q_aprox))
## [1] "El valor apróximado de 'q' es 21"
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
datos <- read.csv("Serie_Examen.csv", header = TRUE, sep = ";")
serie_temporal <- ts(datos$CO2, start = c(1959, 1), frequency = 12)
serie_diferenciada <- diff(serie_temporal, differences = d)
acf_serie_diferenciada <- acf(serie_diferenciada, main = "Función de autocorrelación (ACF)")
pacf_serie_diferenciada <- pacf(serie_diferenciada, main = "Función de autocorrelación parcial (PACF)")
#El código anterior calcula el valor de q aproximado basado en la función de autocorrelación. Esto con una diferencia de 1. El valor aproximado de 1 muestra una estimación basado en p y d.
– Comentario
–
Mostrar pacf, comente el significado de la gráficas
pacf_values <- Pacf(diff(mi_serie_co, differences = d), lag.max = 20)
plot(pacf_values, main = "Función de Autocorrelación Parcial (PACF) de la serie diferenciada")
#El código anterior calcula el valor de q aproximado basado en la función de autocorrelación. Esto con una diferencia de 1. El valor aproximado de 1 muestra una estimación basado en p y d.
– Comentario
–
Verifica con un auto.arima, cuales son los valores de p,d y q, anotalos al final
modelo_auto_arima <- auto.arima(mi_serie_co)
p_valor <- modelo_auto_arima$arma[1]
d_valor <- modelo_auto_arima$arma[6]
q_valor <- modelo_auto_arima$arma[2]
print(paste("El valor de p es:", p_valor))
## [1] "El valor de p es: 1"
print(paste("El valor de d es:", d_valor))
## [1] "El valor de d es: 0"
print(paste("El valor de q es:", q_valor))
## [1] "El valor de q es: 1"
Comentarios
Desarrolla el modelo ARIMA con los datos del punto anterior, carga los resultados en una variable con nombre ‘mi_modelo’
p <- 1
d <- 0
q <- 1
mi_modelo <- arima(mi_serie_co, order = c(p, d, q))
summary(mi_modelo)
##
## Call:
## arima(x = mi_serie_co, order = c(p, d, q))
##
## Coefficients:
## ar1 ma1 intercept
## 0.9675 0.647 339.2790
## s.e. 0.0184 0.046 3.1773
##
## sigma^2 estimated as 0.8689: log likelihood = -219.27, aic = 446.55
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.03896931 0.9321704 0.7818218 0.0106998 0.2304365 0.7029981
## ACF1
## Training set 0.3134469
Revisa sus residuales, comparalos con el proceso estacionario (ruido blanco) que encontraste en el paso 2 y muestra su valor AIC.
residuos <- mi_modelo$residuals^2
plot(mi_modelo$residuals)
plot(residuos)
acf(mi_modelo$residuals)
acf(residuos)
pacf(mi_modelo$residuals)
pacf(residuos)
#La gráfica de los residuos del modelo ARIMA muestra movimientos volátiles dentro de un canal definido por las resistencias y los soportes. En cuanto a la función de autocorrelación, revela que las disminuciones en los residuos son graduales, mientras que las aumentos son más bruscos y definidos. Esta caracterÃstica sugiere una posible correlación en ciertos rezagos, siendo importante explorar más a fondo el significado de este patrón dentro del contexto del modelo ARIMA.
Comentarios
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_serie_co, h = 60)
proyecciones
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jun 1972 345.4803 345.0888 345.8718 344.8816 346.0791
## Jul 1972 346.8467 346.4190 347.2744 346.1925 347.5009
## Aug 1972 348.1831 347.7217 348.6445 347.4775 348.8887
## Sep 1972 349.0450 348.5524 349.5375 348.2917 349.7982
## Oct 1972 349.7977 349.2758 350.3196 348.9995 350.5959
## Nov 1972 350.9214 350.3713 351.4716 350.0801 351.7628
## Dec 1972 352.2955 351.7181 352.8729 351.4125 353.1785
## Jan 1973 352.7984 352.1954 353.4013 351.8763 353.7205
## Feb 1973 352.1693 351.5426 352.7960 351.2108 353.1278
## Mar 1973 350.5279 349.8790 351.1767 349.5355 351.5202
## Apr 1973 348.5115 347.8417 349.1814 347.4871 349.5360
## May 1973 346.9715 346.2812 347.6618 345.9158 348.0272
## Jun 1973 346.8747 346.1490 347.6005 345.7648 347.9846
## Jul 1973 348.2411 347.4948 348.9874 347.0997 349.3825
## Aug 1973 349.5775 348.8110 350.3440 348.4053 350.7498
## Sep 1973 350.4394 349.6534 351.2254 349.2373 351.6415
## Oct 1973 351.1921 350.3871 351.9972 349.9609 352.4234
## Nov 1973 352.3159 351.4919 353.1399 351.0557 353.5761
## Dec 1973 353.6899 352.8471 354.5327 352.4009 354.9789
## Jan 1974 354.1928 353.3320 355.0536 352.8763 355.5093
## Feb 1974 353.5637 352.6858 354.4415 352.2211 354.9063
## Mar 1974 351.9223 351.0284 352.8162 350.5552 353.2894
## Apr 1974 349.9060 348.9966 350.8153 348.5152 351.2967
## May 1974 348.3659 347.4413 349.2906 346.9518 349.7801
## Jun 1974 348.2692 347.3175 349.2208 346.8137 349.7246
## Jul 1974 349.6355 348.6677 350.6033 348.1554 351.1157
## Aug 1974 350.9720 349.9882 351.9557 349.4674 352.4765
## Sep 1974 351.8338 350.8345 352.8331 350.3055 353.3621
## Oct 1974 352.5866 351.5719 353.6012 351.0348 354.1383
## Nov 1974 353.7103 352.6803 354.7403 352.1350 355.2856
## Dec 1974 355.0843 354.0389 356.1298 353.4855 356.6832
## Jan 1975 355.5872 354.5269 356.6475 353.9657 357.2087
## Feb 1975 354.9581 353.8838 356.0324 353.3151 356.6011
## Mar 1975 353.3167 352.2291 354.4043 351.6534 354.9800
## Apr 1975 351.3004 350.2000 352.4008 349.6175 352.9833
## May 1975 349.7603 348.6472 350.8735 348.0579 351.4628
## Jun 1975 349.6636 348.5276 350.7995 347.9263 351.4009
## Jul 1975 351.0300 349.8802 352.1797 349.2715 352.7884
## Aug 1975 352.3664 351.2029 353.5299 350.5870 354.1458
## Sep 1975 353.2282 352.0513 354.4051 351.4283 355.0282
## Oct 1975 353.9810 352.7908 355.1712 352.1607 355.8012
## Nov 1975 355.1047 353.9011 356.3083 353.2640 356.9454
## Dec 1975 356.4788 355.2617 357.6958 354.6174 358.3401
## Jan 1976 356.9816 355.7516 358.2116 355.1005 358.8628
## Feb 1976 356.3525 355.1102 357.5948 354.4526 358.2525
## Mar 1976 354.7111 353.4572 355.9650 352.7935 356.6288
## Apr 1976 352.6948 351.4297 353.9599 350.7600 354.6296
## May 1976 351.1548 349.8784 352.4311 349.2028 353.1068
## Jun 1976 351.0580 349.7615 352.3545 349.0752 353.0408
## Jul 1976 352.4244 351.1156 353.7332 350.4227 354.4260
## Aug 1976 353.7608 352.4397 355.0819 351.7403 355.7813
## Sep 1976 354.6226 353.2895 355.9558 352.5837 356.6616
## Oct 1976 355.3754 354.0303 356.7205 353.3182 357.4326
## Nov 1976 356.4991 355.1419 357.8563 354.4235 358.5748
## Dec 1976 357.8732 356.5038 359.2426 355.7789 359.9675
## Jan 1977 358.3760 356.9949 359.7572 356.2638 360.4883
## Feb 1977 357.7470 356.3547 359.1392 355.6178 359.8762
## Mar 1977 356.1056 354.7029 357.5082 353.9604 358.2507
## Apr 1977 354.0892 352.6765 355.5020 351.9286 356.2499
## May 1977 352.5492 351.1263 353.9721 350.3730 354.7253
plot(proyecciones, main = "Predicción para los siguientes 5 años")