Examen Final

chkPkg <- function(librerias){
  for (lib in librerias) {
    if (!require(lib, character.only = TRUE)) {
      install.packages(lib)
      if (!require(lib, character.only = TRUE)) {
        stop(paste("No se pudo cargar la librería:", lib))
      }
    }
    library(lib, character.only = TRUE)
  }
}

library("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.5.1     ✔ 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

library("forecast")

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library("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

library("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

library("MASS")

## 
## Attaching package: 'MASS'
## 
## The following object is masked from 'package:dplyr':
## 
##     select

library("tseries")
library("stats")
library("car")
library("lmtest")
library("urca")
library("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

library("openxlsx")

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 Temporal CO2", ylab = "CO2", xlab = "Tiempo", col = "blue")

#Tendencia: Existe una clara tendencia creciente en la concentración de CO2 a través del tiempo, mostrando un aumento sostenido desde 1959 hasta aproximadamente 1972.

#Estacionalidad: Se observa un patrón cíclico, repetitivo cada año, indicando estacionalidad mensual muy marcada, con picos y valles que probablemente responden a factores estacionales anuales.

#Estacionariedad: Debido a la presencia de una tendencia creciente evidente, esta serie no es estacionaria.

acf(mi_serie_co, main = "Función de Autocorrelación (ACF) - CO2")

#La gráfica ACF muestra una serie temporal con una fuerte dependencia temporal, tendencia creciente y estacionalidad clara. Por tanto, podemos concluir que la serie de CO2 no es estacionaria en su estado actual.Por lo que es necesario realizar una diferenciación.

descomposicion_stl <- stl(mi_serie_co, s.window = "periodic")
plot(descomposicion_stl, main = "Descomposición STL de la Serie CO2")

# Tendencia: Se observa una clara tendencia creciente en el componente “trend” (tendencia). Esto confirma visualmente lo indicado anteriormente, que la serie presenta un aumento progresivo en los niveles de CO2 con el tiempo.

Estacionalidad: El componente “seasonal” muestra un patrón cíclico claro y regular, repetido anualmente, confirmando la fuerte estacionalidad anual observada previamente tanto visualmente como en el análisis ACF.

Residuales (remainder): El componente residual muestra un comportamiento aleatorio alrededor de cero, sin patrones claros, indicando que los componentes de tendencia y estacionalidad capturan adecuadamente la estructura de la serie.

Estacionariedad:Dado que hay una tendencia marcada y una estacionalidad definida, se confirma nuevamente que la serie original no es estacionaria.

#La descomposición STL respalda plenamente las conclusiones anteriores.

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

# Limpiar información para evitar errores
mi_serie_co <- na.omit(mi_serie_co)
ndiffs(mi_serie_co)

## [1] 1

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

Comentarios: La serie original no es estacionaria “pura”, pero es estacionaria alrededor de una tendencia, por lo que es necesario realizar una diferenciación (d=1) para hacerla completamente estacionaria y proceder con el modelado ARIMA.

plot(diff(mi_serie_co, differences=1), main="Primera Diferencia")

#Comentarios: La primera diferencia elimina exitosamente la tendencia pero no la estacionalidad. La serie ahora es más estacionaria, pero aún requiere considerar la estacionalidad.

adf.test(diff(mi_serie_co, differences=1))

## 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.01

#Se confirma claramente que el valor de diferenciación adecuado es d = 1.

d <- ndiffs(mi_serie_co)

mi_serie_cod1 <- diff(mi_serie_co, difference=d)

MRL <- lm(mi_serie_cod1 ~time(mi_serie_cod1))
summary(MRL)

## 
## Call:
## lm(formula = mi_serie_cod1 ~ time(mi_serie_cod1))
## 
## 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(mi_serie_cod1)  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.9486

El modelo lineal muestra una pendiente (coeficiente) no significativa (p=0.949), indicando que, después de la diferenciación, ya no existe tendencia lineal significativa en la serie. La diferenciación fue efectiva, dejando residuos sin patrón claro (media cerca de cero y errores distribuidos).

Identifica los valor p y q

pacf(mi_serie_cod1)

# P=1

acf(mi_serie_cod1)

Q=1

p <- 1
d <- 1
q <- 1

modelo.arima1 <- arima(mi_serie_co, order=c(p,d,q))

coef(modelo.arima1)

##       ar1       ma1 
## 0.5832764 0.3376746

plot(modelo.arima1$residuals)

acf(modelo.arima1$residuals)

pacf(modelo.arima1$residuals)

modelo.arima1$aic

## [1] 402.6446

par(mfrow = c(1,2))
acf(mi_serie_co)
pacf(mi_serie_co)

#Comentarios: La gráfica de autocorrelación (ACF) muestra una caída abrupta después del primer rezago, sugiriendo claramente un modelo de media móvil con parámetro q=1. Por otro lado, la gráfica de autocorrelación parcial (PACF) presenta un corte significativo inmediatamente después del primer rezago, indicando un modelo autorregresivo con parámetro p=1. En resumen, ambas gráficas sugieren de manera conjunta el uso de un modelo ARIMA(1,1,1).

auto.arima(mi_serie_co, ic="aic", trace=T)

## 
##  Fitting models using approximations to speed things up...
## 
##  ARIMA(2,0,2)(1,1,1)[12] with drift         : 77.04773
##  ARIMA(0,0,0)(0,1,0)[12] with drift         : 211.3175
##  ARIMA(1,0,0)(1,1,0)[12] with drift         : 100.3196
##  ARIMA(0,0,1)(0,1,1)[12] with drift         : 136.3408
##  ARIMA(0,0,0)(0,1,0)[12]                    : 534.7961
##  ARIMA(2,0,2)(0,1,1)[12] with drift         : 73.8594
##  ARIMA(2,0,2)(0,1,0)[12] with drift         : 115.8998
##  ARIMA(2,0,2)(0,1,2)[12] with drift         : 75.02558
##  ARIMA(2,0,2)(1,1,0)[12] with drift         : 88.54737
##  ARIMA(2,0,2)(1,1,2)[12] with drift         : 71.14947
##  ARIMA(2,0,2)(2,1,2)[12] with drift         : 62.27961
##  ARIMA(2,0,2)(2,1,1)[12] with drift         : 67.28559
##  ARIMA(1,0,2)(2,1,2)[12] with drift         : 58.73856
##  ARIMA(1,0,2)(1,1,2)[12] with drift         : 70.05191
##  ARIMA(1,0,2)(2,1,1)[12] with drift         : 64.08053
##  ARIMA(1,0,2)(1,1,1)[12] with drift         : 74.82702
##  ARIMA(0,0,2)(2,1,2)[12] with drift         : 86.74447
##  ARIMA(1,0,1)(2,1,2)[12] with drift         : 57.14473
##  ARIMA(1,0,1)(1,1,2)[12] with drift         : 69.11954
##  ARIMA(1,0,1)(2,1,1)[12] with drift         : 63.56778
##  ARIMA(1,0,1)(1,1,1)[12] with drift         : 74.8631
##  ARIMA(0,0,1)(2,1,2)[12] with drift         : 108.1426
##  ARIMA(1,0,0)(2,1,2)[12] with drift         : 73.04009
##  ARIMA(2,0,1)(2,1,2)[12] with drift         : 60.32313
##  ARIMA(0,0,0)(2,1,2)[12] with drift         : 142.667
##  ARIMA(2,0,0)(2,1,2)[12] with drift         : 60.00372
##  ARIMA(1,0,1)(2,1,2)[12]                    : Inf
## 
##  Now re-fitting the best model(s) without approximations...
## 
##  ARIMA(1,0,1)(2,1,2)[12] with drift         : 84.92471
## 
##  Best model: ARIMA(1,0,1)(2,1,2)[12] with drift

## 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.96

arima.opt <- arima(mi_serie_co, order=c(1,0,1), seasonal = list(order = c(2, 1, 2), period = 12))
summary(arima.opt)

## 
## Call:
## arima(x = mi_serie_co, order = c(1, 0, 1), seasonal = list(order = c(2, 1, 2), 
##     period = 12))
## 
## Coefficients:
##          ar1      ma1     sar1     sar2     sma1     sma2
##       0.9997  -0.5023  -0.5138  -0.0569  -0.1965  -0.4410
## s.e.  0.0006   0.0735   0.3272   0.1549   0.3282   0.2878
## 
## sigma^2 estimated as 0.09081:  log likelihood = -39.83,  aic = 93.65
## 
## Training set error measures:
##                      ME      RMSE       MAE       MPE       MAPE      MASE
## Training set 0.04920222 0.3035994 0.2498824 0.0147126 0.07362382 0.2246891
##                    ACF1
## Training set 0.05206132

checkresiduals(arima.opt)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,1)(2,1,2)[12]
## Q* = 18.457, df = 18, p-value = 0.426
## 
## Model df: 6.   Total lags used: 24

coef(arima.opt)

##         ar1         ma1        sar1        sar2        sma1        sma2 
##  0.99968618 -0.50226129 -0.51377818 -0.05694573 -0.19650821 -0.44102381

arima.opt$aic

## [1] 93.65198

#El mejor modelo arima es ARIMA(1,0,1)(2,1,2)[12] con un AIC del 93.65198

preds <- predict(arima.opt, n.ahead = 5*12)
ts.plot(cbind(mi_serie_co, preds$pred), lty=c(1,3))

preds$pred

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 1972                                              345.4190 346.8200 348.1627
## 1973 352.9395 352.3244 350.6199 348.4795 347.0048 346.9089 348.2975 349.6737
## 1974 354.3092 353.7069 352.0046 349.9240 348.3767 348.2547 349.6500 350.9990
## 1975 355.7414 355.1230 353.4111 351.2884 349.8024 349.6809 351.0727 352.4331
## 1976 357.1252 356.5136 354.8059 352.7007 351.1867 351.0657 352.4582 353.8136
## 1977 358.5221 357.9072 356.1972 354.0847 352.5809                           
##           Sep      Oct      Nov      Dec
## 1972 349.0322 349.7497 350.9511 352.2753
## 1973 350.4974 351.2209 352.2630 353.7178
## 1974 351.8563 352.5631 353.7017 355.1056
## 1975 353.2751 353.9894 355.0868 356.5087
## 1976 354.6609 355.3716 356.4840 357.8989
## 1977