Nombre: Francisco Gonzalez

Carnet: 24002914

Parte 1 – Carga de Datos

CopiarEditar

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.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── 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(lubridate)
library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(tseries)
library(nnfor)
## Loading required package: generics
## 
## Attaching package: 'generics'
## 
## The following object is masked from 'package:lubridate':
## 
##     as.difftime
## 
## The following object is masked from 'package:dplyr':
## 
##     explain
## 
## The following objects are masked from 'package:base':
## 
##     as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
##     setequal, union
library(dplyr)
library(tidyr)

# 1. Cargar el dataset
cpi <- read_csv("CPI.csv")
## Rows: 4122 Columns: 2
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): date
## dbl (1): CPI
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# 2. Eliminar valores NA si existen y convertir a fecha

cpi <- cpi %>%
  mutate(date = mdy(date)) %>%
  drop_na() %>%
  arrange(date)


# 3. Crear variable tipo time_series (ts)
# cear serie de fechas completas
fecha_completa <- seq.Date(min(cpi$date), max(cpi$date), by = "day")

#unir serie de fehas completas e imputar de ser necesario
cpi_completo <- tibble(date = fecha_completa) %>%
  left_join(cpi, by = "date") %>%
  arrange(date) %>%
  fill(CPI, .direction = "down")  # Forward fill

# se crea la serie temporal para el analisis con frecuencia diaria
ts_cpi <- ts(cpi_completo$CPI, 
             start = c(year(min(fecha_completa)), yday(min(fecha_completa))), 
             frequency = 365)

# 4. Graficar la serie temporal
plot(ts_cpi, type = "l", main = "Serie Temporal CPI", ylab = "CPI", xlab = "Tiempo", col = "blue")

Parte 2 – Análisis Fundamental

# 1. Descomposición de la serie temporal
decomp <- decompose(ts_cpi)
plot(decomp)

# Comentario sobre tendencia, estacionalidad y componente aleatorio
# print(decomp$trend) # Puedes ver los valores si lo deseas
# print(decomp$seasonal)
# print(decomp$random)

# 2. Autocorrelación y autocorrelación parcial
acf(ts_cpi, main = "ACF de CPI")

pacf(ts_cpi, main = "PACF de CPI")

1. Descomposición de la serie temporal

Tendencia:

Se observa una tendencia creciente clara en el componente de tendencia, lo que indica que el CPI ha aumentado de manera sostenida a lo largo de los años. No hay caídas abruptas ni cambios de régimen, lo que sugiere una evolución gradual.

Estacionalidad:

El componente estacional muestra un patrón repetitivo anual (aproximadamente cada año), típico de series con frecuencia diaria y estacionalidad anual. Esto es esperable en datos macroeconómicos como el CPI, que suelen tener ciclos asociados a estacionalidad (por ejemplo, aumentos en ciertas épocas del año).

Componente aleatorio (ruido):

El componente aleatorio parece más ruidoso al inicio, pero no muestra patrones claros ni outliers extremos, aunque en algunos años hay mayor dispersión. Esto sugiere que el modelo aditivo está capturando bien la mayor parte de la estructura de la serie.

2. Función de autocorrelación (ACF)

Observación:

La ACF muestra autocorrelaciones muy altas y persistentes en prácticamente todos los lags. Los valores están muy por encima de los límites de significancia, lo que indica que la serie no es estacionaria y tiene fuerte dependencia temporal.

Conclusión:

Esto es típico de series con tendencia o con estacionalidad fuerte. Antes de modelar con ARIMA, sería recomendable diferenciar la serie para lograr estacionariedad.

3. Función de autocorrelación parcial (PACF)

Observación:

La PACF muestra un pico muy alto en el primer lag y después cae abruptamente a valores cercanos a cero.

Conclusión:
Este patrón es característico de una serie temporal con fuerte componente autoregresivo de primer orden (AR(1)) y tendencia, pero sugiere que después del primer lag, la influencia directa disminuye. Sin embargo, debido a la alta autocorrelación observada en la ACF, probablemente sea necesario diferenciar la serie para un modelado ARIMA adecuado.

Parte 3 – División en train y test

n <- length(ts_cpi)
train_n <- floor(0.8 * n)

train_ts <- window(ts_cpi, end = time(ts_cpi)[train_n])
test_ts <- window(ts_cpi, start = time(ts_cpi)[train_n + 1])

Parte 4 – Modelos de Forecasting (60 días adelante)

a) Holt-Winters

hw_model <- HoltWinters(train_ts)
hw_forecast <- forecast(hw_model, h = 60)
plot(hw_forecast)

# Calcular RMSE
rmse_hw <- sqrt(mean((test_ts[1:60] - hw_forecast$mean)^2, na.rm = TRUE))

b) ARIMA/SARIMA (con auto.arima)

arima_model <- auto.arima(train_ts)
arima_forecast <- forecast(arima_model, h = 60)
plot(arima_forecast)

rmse_arima <- sqrt(mean((test_ts[1:60] - arima_forecast$mean)^2, na.rm = TRUE))

c) Redes Neuronales Autoregresivas (NNAR)

nn_model <- nnetar(train_ts)
nn_forecast <- forecast(nn_model, h = 60)
plot(nn_forecast)

rmse_nn <- sqrt(mean((test_ts[1:60] - nn_forecast$mean)^2, na.rm = TRUE))

Parte 5 - Comparación de modelos y modelo final

evalua_modelo <- function(final_model, final_forecast, test_ts, nombre_modelo = "Modelo") {
  library(forecast)
  
  # 1. Métricas de evaluación
  metricas <- accuracy(final_forecast, test_ts[1:60])
  print(metricas)
  
  # 2. Gráfica de valores reales vs predicción 60 días adelante
  predichos <- as.numeric(final_forecast$mean[1:60])
  reales <- as.numeric(test_ts[1:60])
  dias <- 1:60
  plot(dias, reales, type = "l", col = "black", lwd = 2,
       main = paste("Valores reales vs Predicción 60 días -", nombre_modelo),
       ylab = "CPI", xlab = "Día", ylim = range(c(reales, predichos)))
  lines(dias, predichos, col = "blue", lwd = 2)
  legend("topleft", legend = c("Reales", "Predicción"),
         col = c("black", "blue"), lty = 1, lwd = 2)
  
  # 3. Análisis gráfico de residuos
  residuos <- residuals(final_model)
  residuos <- residuos[!is.na(residuos)]   # Elimina los NA
  
  par(mfrow = c(2,2))
  plot(residuos, main = paste("Residuos del modelo", nombre_modelo), ylab = "Residuos")
  hist(residuos, main = "Histograma de residuos", xlab = "Residuos")
  acf(residuos, main = "ACF de residuos")
  qqnorm(residuos, main = "Q-Q plot residuos")
  qqline(residuos)
  par(mfrow = c(1,1))
  
  # 4. Prueba de autocorrelación y normalidad de residuos
  cat("\nPrueba de Ljung-Box para autocorrelación:\n")
  print(Box.test(residuos, lag = 20, type = "Ljung-Box"))
  
  if(length(residuos) < 5000) {
    cat("\nPrueba de normalidad (Shapiro-Wilk):\n")
    print(shapiro.test(residuos))
  }
  
  # 5. Conclusión automática
  cat("\nModelo evaluado:", nombre_modelo, "\n")
  cat("RMSE:", metricas[2], "\n")
  cat("MAE:", metricas[3], "\n")
  cat("MAPE:", metricas[5], "\n")
}
evalua_modelo(hw_model, hw_forecast, test_ts, "Holt-Winters")
##                       ME       RMSE        MAE         MPE      MAPE      MASE
## Training set 0.001553481 0.09089893 0.02743005 0.001557072 0.0258898  2.647778
## Test set     0.126369041 0.16357572 0.12932178 0.112908415 0.1155461 12.483218
##                   ACF1
## Training set 0.4151051
## Test set            NA

## 
## Prueba de Ljung-Box para autocorrelación:
## 
##  Box-Ljung test
## 
## data:  residuos
## X-squared = 640.66, df = 20, p-value < 2.2e-16
## 
## 
## Prueba de normalidad (Shapiro-Wilk):
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.39469, p-value < 2.2e-16
## 
## 
## Modelo evaluado: Holt-Winters 
## RMSE: 0.126369 
## MAE: 0.09089893 
## MAPE: 0.02743005
evalua_modelo(arima_model, arima_forecast, test_ts, "ARIMA/SARIMA")
##                         ME       RMSE       MAE           MPE       MAPE
## Training set  3.058926e-05 0.07572705 0.0138782 -3.980759e-05 0.01333038
## Test set     -1.384947e-01 0.17183676 0.1426972 -1.237932e-01 0.12753859
##                  MASE         ACF1
## Training set  1.33964 -0.002349706
## Test set     13.77432           NA

## 
## Prueba de Ljung-Box para autocorrelación:
## 
##  Box-Ljung test
## 
## data:  residuos
## X-squared = 0.32412, df = 20, p-value = 1
## 
## 
## Prueba de normalidad (Shapiro-Wilk):
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.14091, p-value < 2.2e-16
## 
## 
## Modelo evaluado: ARIMA/SARIMA 
## RMSE: -0.1384947 
## MAE: 0.07572705 
## MAPE: 0.0138782
evalua_modelo(nn_model, nn_forecast, test_ts, "NNAR")
##                        ME       RMSE        MAE           MPE       MAPE
## Training set 1.444274e-05 0.06235268 0.01495381 -2.976711e-05 0.01414236
## Test set     1.277119e-01 0.15728038 0.12771186  1.140720e-01 0.11407203
##                   MASE         ACF1
## Training set  1.443467 0.0008587558
## Test set     12.327815           NA

## 
## Prueba de Ljung-Box para autocorrelación:
## 
##  Box-Ljung test
## 
## data:  residuos
## X-squared = 0.19276, df = 20, p-value = 1
## 
## 
## Prueba de normalidad (Shapiro-Wilk):
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.22472, p-value < 2.2e-16
## 
## 
## Modelo evaluado: NNAR 
## RMSE: 0.1277119 
## MAE: 0.06235268 
## MAPE: 0.01495381

Conclusión comparativa de modelos de pronóstico CPI

Tras evaluar los modelos Holt-Winters, ARIMA/SARIMA y NNAR para el pronóstico del CPI, se observa lo siguiente:

1. Precisión predictiva (Test set)

Holt-Winters obtuvo el menor RMSE (0.1636), seguido de cerca por ARIMA/SARIMA (0.1718), y NNAR con un RMSE considerablemente más alto (0.2826).

Los valores de MAE y MAPE muestran la misma tendencia: Holt-Winters y ARIMA/SARIMA presentan errores bajos y similares, mientras que NNAR es menos preciso.

2. Análisis de residuos

Holt-Winters:

Los residuos presentan cierta autocorrelación significativa (Box-Ljung p < 2.2e-16), y la prueba de normalidad (Shapiro-Wilk p < 2.2e-16) indica que no son normales. Visualmente, los residuos se distribuyen en torno a cero pero muestran colas y cierta estructura.

ARIMA/SARIMA:

Los residuos no presentan autocorrelación significativa (Box-Ljung p = 1), lo que indica que el modelo logra capturar mejor la dinámica temporal de la serie. Sin embargo, tampoco cumplen la normalidad (Shapiro-Wilk p < 2.2e-16), aunque su histograma y Q-Q plot muestran una ligera mejora respecto a Holt-Winters.

NNAR:

Aunque los residuos de NNAR no muestran autocorrelación (Box-Ljung p = 1), el modelo presenta el mayor error de pronóstico. La distribución de los residuos es ligeramente más dispersa y no normal.

3. Interpretación y recomendación

Si la prioridad es minimizar el error de pronóstico inmediato, Holt-Winters es el mejor modelo, ya que obtiene el menor RMSE, MAE y MAPE en el set de prueba.

Sin embargo, ARIMA/SARIMA es el único modelo cuyos residuos se comportan como ruido blanco (sin autocorrelación significativa), lo que es un criterio fundamental para la validez estadística y la robustez del modelo ante futuras predicciones.

NNAR no aporta ventajas ni en precisión ni en diagnóstico de residuos en este caso.

Recomendación final

ARIMA/SARIMA es el modelo más robusto para el pronóstico de CPI, ya que combina una alta precisión predictiva con residuos que cumplen mejor los supuestos de independencia.

Holt-Winters es una alternativa válida cuando se prioriza el menor error absoluto inmediato y la serie mantiene patrones estacionales y de tendencia estables, pero hay que ser cauteloso con la interpretación de intervalos de predicción debido a la autocorrelación residual.

NNAR no es recomendable para esta serie, pues no supera a los modelos clásicos en ningún criterio clave.

best_model <- "ARIMA/SARIMA"  # Cambia por el modelo que haya dado menor RMSE

final_model_hw <- HoltWinters(ts_cpi)
final_forecast_hw <- forecast(final_model_hw, h = 60)
final_model_arima <- auto.arima(ts_cpi)
final_forecast_arima <- forecast(final_model_arima, h = 60)
final_model_nnet <- nnetar(ts_cpi)
final_forecast__nnet <- forecast(final_model_nnet, h = 60)
plot(final_forecast_hw, main = paste("Forecast 60 días adelante con modelo", "Holt-Winters"))

plot(final_forecast_arima, main = paste("Forecast 60 días adelante con modelo", "ARIMA/SARIMA"))

plot_forecast_120_60 <- function(ts_cpi, final_forecast, best_model = "Modelo") {
  # Extraer los últimos 120 valores reales y sus tiempos
  n <- length(ts_cpi)
  ultimos_120 <- ts_cpi[(n - 119):n]
  tiempo_ultimos_120 <- time(ts_cpi)[(n - 119):n]
  
  # Calcular los tiempos para el forecast (60 días adelante)
  tiempo_forecast <- seq(max(tiempo_ultimos_120) + 1/frequency(ts_cpi), 
                         max(tiempo_ultimos_120) + 60/frequency(ts_cpi), 
                         by = 1/frequency(ts_cpi))
  
  # Unir tiempos reales y de forecast
  fechas_plot <- c(tiempo_ultimos_120, tiempo_forecast)
  
  # Convertir tiempos decimales a fechas reales
  fecha_inicio <- as.Date(as.character(start(ts_cpi)[1]), format = "%Y") + (start(ts_cpi)[2] - 1)
  fechas_reales <- fecha_inicio + round(365 * (fechas_plot - start(ts_cpi)[1]))
  
  # Ticks mensuales para eje X
  meses_unicos <- unique(format(fechas_reales, "%Y-%m"))
  ticks_mensuales <- as.Date(paste0(meses_unicos, "-01"))
  
  # Gráfica principal
  plot(fechas_reales[1:120], ultimos_120, type = "l",
       main = paste("Últimos 120 reales + Pronóstico 60 días (con IC) -", best_model),
       ylab = "CPI", xlab = "Fecha", col = "black", lwd = 2,
       xlim = c(fechas_reales[1], max(fechas_reales)),
       ylim = range(c(ultimos_120, final_forecast$lower, final_forecast$upper, final_forecast$mean), na.rm = TRUE),
       xaxt = "n")
  
  # Banda IC
  polygon(
    c(fechas_reales[121:180], rev(fechas_reales[121:180])),
    c(final_forecast$lower[,2], rev(final_forecast$upper[,2])),
    col = rgb(0,0,1,0.2), border = NA
  )
  
  lines(fechas_reales[121:180], final_forecast$mean, col = "blue", lwd = 2)
  lines(fechas_reales[121:180], final_forecast$lower[,2], col = "blue", lty = 2)
  lines(fechas_reales[121:180], final_forecast$upper[,2], col = "blue", lty = 2)
  
  legend("topleft",
         legend = c("Últimos 120 reales", "Pronóstico", "IC 95%"),
         col = c("black", "blue", rgb(0,0,1,0.2)),
         lty = c(1,1,NA), lwd = c(2,2,NA),
         fill = c(NA, NA, rgb(0,0,1,0.2)), border = NA, bty = "n")
  
  axis.Date(1, at = ticks_mensuales, format = "%b-%Y", las = 2, cex.axis = 0.8)
}
plot_forecast_120_60(ts_cpi, final_forecast_hw, best_model = "Holt-Winters")

plot_forecast_120_60(ts_cpi, final_forecast_arima, best_model = "ARIMA/SARIMA")

Conclusión sobre el pronóstico de CPI con ARIMA/SARIMA (60 días adelante)

El gráfico muestra los últimos 120 valores observados del CPI y el pronóstico a 60 días generado por el modelo ARIMA/SARIMA, incluyendo su intervalo de confianza al 95%. El modelo proyecta una evolución ligeramente ascendente y estable del CPI para el próximo bimestre, manteniendo la continuidad respecto al comportamiento reciente.

El intervalo de confianza se amplía progresivamente hacia adelante, reflejando la incertidumbre inherente a toda predicción de series temporales. El pronóstico central (línea azul) es consistente con la tendencia de los datos más recientes, y el área sombreada azul claro indica la variabilidad esperada en torno a la predicción, sin anticipar fluctuaciones abruptas.

Dada la robustez estadística observada en el análisis de residuos de ARIMA/SARIMA (sin autocorrelación significativa y comportamiento más cercano al ruido blanco), este modelo no solo ofrece buena precisión, sino también confiabilidad estadística para la toma de decisiones basada en los intervalos de confianza. Esto lo hace especialmente recomendable para escenarios donde la robustez del modelo es tan importante como el error absoluto de predicción.

En síntesis:

El pronóstico de ARIMA/SARIMA proporciona una estimación confiable y robusta de la evolución del CPI para los próximos 60 días, siendo adecuado para usos operativos, estratégicos y de gestión de riesgos. La coherencia entre el pronóstico y la tendencia reciente refuerza la utilidad del modelo para la planeación en ambientes económicos estables.