library(forecast)
## Warning: package 'forecast' was built under R version 4.3.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(tseries)
## Warning: package 'tseries' was built under R version 4.3.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.3.3
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(lubridate)
## Warning: package 'lubridate' was built under R version 4.3.3
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(zoo)
## Warning: package 'zoo' was built under R version 4.3.3
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(xts)
## Warning: package 'xts' was built under R version 4.3.3
## 
## ######################### 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
#Primero cargaremos los datos a través del archivo de Excel y transformar la columna Mes en formato Date

library(readxl)
Base_de_Datos_18_23_Calificadora_de_Riesgo_y_Rentabilidad_ <- read_excel("Base de Datos 18-23 (Calificadora de Riesgo y Rentabilidad).xlsx")

Datos <- Base_de_Datos_18_23_Calificadora_de_Riesgo_y_Rentabilidad_
Datos$Mes <- as.Date(Datos$Mes, format = "%m/%d/%Y")
#Debido a que los datos en el Excel están en porcentaje, escalaremos los valores multiplicandolos por 100

datos_porcentajes <- Datos %>%
  mutate_at(vars(-Mes, -AFP), function(x) as.numeric(x) * 100)

print(datos_porcentajes)
## # A tibble: 504 × 10
##    Mes        AFP       Rentabilidad `C-1`  `C-2` `C-3`   AAA    AA      A   BBB
##    <date>     <chr>            <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl>
##  1 2023-12-01 ROMANA            7.31 18.7  0      0      80.8  0     0.529  0   
##  2 2023-12-01 RESERVAS          8.94  4.33 0.0492 0      72.2  8.27  4.77  10.4 
##  3 2023-12-01 POPULAR           9.58  2.97 0.685  0      71.5  8.45  4.53  11.8 
##  4 2023-12-01 SIEMBRA           8.32  2.90 0.545  0      76.2  9.74  4.58   6.07
##  5 2023-12-01 CRECER            7.81  1.90 0      0      70.0 18.3   2.76   7.06
##  6 2023-12-01 JMMB- BDI         9.94  1.29 5.29   0      71.6  1.60 16.3    3.95
##  7 2023-12-01 ATLÁNTICO         7.72  0    0.987  0.570  76.6  2.90 12.9    6.04
##  8 2023-11-01 ROMANA            7.35 19.3  0      0      80.2  0     0.528  0   
##  9 2023-11-01 RESERVAS          9.05  3.12 0.0715 0      75.9  6.95  4.63   9.32
## 10 2023-11-01 SIEMBRA           8.95  2.90 0.505  0      77.2  8.85  4.63   5.94
## # ℹ 494 more rows
#Para la trasformación de los datos a series temporales, debemos separar los fondos del data set principal de datos, unicamente trabajaremos con la Rentabilidad, por lo que separaremos de cada fondo el Mes y la Rentabilidad.

datos_atlantico <- subset(datos_porcentajes, AFP == "ATLÁNTICO", select = c(Mes, Rentabilidad))
datos_crecer <- subset(datos_porcentajes, AFP == "CRECER", select = c(Mes, Rentabilidad))
datos_Jmmb <- subset(datos_porcentajes, AFP == "JMMB- BDI", select = c(Mes, Rentabilidad))
datos_popular <- subset(datos_porcentajes, AFP == "POPULAR", select = c(Mes, Rentabilidad))
datos_reservas <- subset(datos_porcentajes, AFP == "RESERVAS", select = c(Mes, Rentabilidad))
datos_romana <- subset(datos_porcentajes, AFP == "ROMANA", select = c(Mes, Rentabilidad))
datos_siembra <- subset(datos_porcentajes, AFP == "SIEMBRA", select = c(Mes, Rentabilidad))
#Para crear las series temporales, se crearan columna por columna las series y luego se combinarán

# Obtener la fecha mínima y máxima de la columna 'Mes'
fecha_minima <- min(Datos$Mes)
fecha_maxima <- max(Datos$Mes)

# Crear las series temporales para cada columna
ATLANTICO <- ts(datos_atlantico$Rentabilidad, start = c(year(fecha_minima), month(fecha_minima)), end = c(year(fecha_maxima), month(fecha_maxima)), frequency = 12)

SIEMBRA <- ts(datos_siembra$Rentabilidad, start = c(year(fecha_minima), month(fecha_minima)), end = c(year(fecha_maxima), month(fecha_maxima)), frequency = 12)

CRECER <- ts(datos_crecer$Rentabilidad, start = c(year(fecha_minima), month(fecha_minima)), end = c(year(fecha_maxima), month(fecha_maxima)), frequency = 12)

ROMANA <- ts(datos_romana$Rentabilidad, start = c(year(fecha_minima), month(fecha_minima)), end = c(year(fecha_maxima), month(fecha_maxima)), frequency = 12)

RESERVAS <- ts(datos_reservas$Rentabilidad, start = c(year(fecha_minima), month(fecha_minima)), end = c(year(fecha_maxima), month(fecha_maxima)), frequency = 12)

POPULAR <- ts(datos_popular$Rentabilidad, start = c(year(fecha_minima), month(fecha_minima)), end = c(year(fecha_maxima), month(fecha_maxima)), frequency = 12)

JMMB_BDI <- ts(datos_Jmmb$Rentabilidad, start = c(year(fecha_minima), month(fecha_minima)), end = c(year(fecha_maxima), month(fecha_maxima)), frequency = 12)
#Combinando las series en un mismo Data set, para fines de graficar las series temporales

Datos_Serie_temporal <- cbind(ATLANTICO,CRECER,JMMB_BDI,POPULAR,RESERVAS,ROMANA,SIEMBRA)
plot(Datos_Serie_temporal)

#Grafica unificada

plot(datos_atlantico$Mes,datos_atlantico$Rentabilidad, "l", col = "Blue",
     main = "Rentabilidad Atlantico", xlab = "Fecha", ylab = "Rentabilidad")
lines(datos_crecer$Mes, datos_crecer$Rentabilidad,"l", col = "red")
lines(datos_Jmmb$Mes, datos_Jmmb$Rentabilidad,"l", col = "green")
lines(datos_popular$Mes, datos_popular$Rentabilidad,"l", col = "gray")
lines(datos_reservas$Mes, datos_reservas$Rentabilidad,"l", col = "orange")
lines(datos_romana$Mes, datos_romana$Rentabilidad,"l", col = "Black")
lines(datos_siembra$Mes, datos_siembra$Rentabilidad,"l", col = "purple")

  legend("topright", legend = c("Atlantico", "Crecer", "Jmmb - BDI", "Popular", "Reservas", "Romana", "Siembre"), col = c("blue", "red", "green", "gray","orange","Black","purple"), lty = 1)

#Ahora estadisticamente veremos si estas series son estacionarias o no estacionarias con la prueba Dickey-Fuller

apply_adf_tests <- function(series_list, series_names) {
  if (length(series_list) != length(series_names)) {
    stop("La longitud de series_list y series_names debe ser la misma")
  }

  results <- list()

  for (i in 1:length(series_list)) {
    series <- series_list[[i]]
    name <- series_names[i]
    test_result <- adf.test(series)
    results[[name]] <- test_result
    cat("Prueba Dickey-Fuller para", name, "\n")
    print(test_result)
    cat("\n")
  }

  return(results)
}
series_list <- list(ATLANTICO,CRECER,SIEMBRA,ROMANA,RESERVAS,POPULAR,JMMB_BDI)
series_names <- c("ATLANTICO", "CRECER","SIEMBRA","ROMANA","RESERVAS","POPULAR","JMMB_BDI")
Prueba_Dickey_Fuller <- apply_adf_tests(series_list, series_names)
## Prueba Dickey-Fuller para ATLANTICO 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  series
## Dickey-Fuller = -2.1987, Lag order = 4, p-value = 0.4944
## alternative hypothesis: stationary
## 
## 
## Prueba Dickey-Fuller para CRECER 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  series
## Dickey-Fuller = -1.8888, Lag order = 4, p-value = 0.6206
## alternative hypothesis: stationary
## 
## 
## Prueba Dickey-Fuller para SIEMBRA 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  series
## Dickey-Fuller = -2.642, Lag order = 4, p-value = 0.3139
## alternative hypothesis: stationary
## 
## 
## Prueba Dickey-Fuller para ROMANA 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  series
## Dickey-Fuller = -3.549, Lag order = 4, p-value = 0.044
## alternative hypothesis: stationary
## 
## 
## Prueba Dickey-Fuller para RESERVAS 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  series
## Dickey-Fuller = -2.7278, Lag order = 4, p-value = 0.279
## alternative hypothesis: stationary
## 
## 
## Prueba Dickey-Fuller para POPULAR 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  series
## Dickey-Fuller = -2.5428, Lag order = 4, p-value = 0.3543
## alternative hypothesis: stationary
## 
## 
## Prueba Dickey-Fuller para JMMB_BDI 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  series
## Dickey-Fuller = -2.5589, Lag order = 4, p-value = 0.3478
## alternative hypothesis: stationary
apply_adf_tests2 <- function(series_list, series_names) {
  if (length(series_list) != length(series_names)) {
    stop("La longitud de series_list y series_names debe ser la misma")
  }

  results <- list()

  for (i in 1:length(series_list)) {
    series <- series_list[[i]]
    name <- series_names[i]
    test_result <- adf.test(series)
    results[[name]] <- test_result
    cat("Prueba Dickey-Fuller para", name, "\n")
    print(test_result)
    if (test_result$p.value < 0.05) {
      cat(name, "es estacionaria (se rechaza la hipótesis nula)\n")
    } else {
      cat(name, "es no estacionaria (no se rechaza la hipótesis nula)\n")
    }
    cat("\n")
  }

  return(results)
}
#Automaticemos el proceso de rechazar la hipotesis nula nos diga si es estacionaria y en caso contrario decir que es estacionaria

Prueba_Dickey_Fuller2 <- apply_adf_tests2(series_list, series_names)
## Prueba Dickey-Fuller para ATLANTICO 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  series
## Dickey-Fuller = -2.1987, Lag order = 4, p-value = 0.4944
## alternative hypothesis: stationary
## 
## ATLANTICO es no estacionaria (no se rechaza la hipótesis nula)
## 
## Prueba Dickey-Fuller para CRECER 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  series
## Dickey-Fuller = -1.8888, Lag order = 4, p-value = 0.6206
## alternative hypothesis: stationary
## 
## CRECER es no estacionaria (no se rechaza la hipótesis nula)
## 
## Prueba Dickey-Fuller para SIEMBRA 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  series
## Dickey-Fuller = -2.642, Lag order = 4, p-value = 0.3139
## alternative hypothesis: stationary
## 
## SIEMBRA es no estacionaria (no se rechaza la hipótesis nula)
## 
## Prueba Dickey-Fuller para ROMANA 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  series
## Dickey-Fuller = -3.549, Lag order = 4, p-value = 0.044
## alternative hypothesis: stationary
## 
## ROMANA es estacionaria (se rechaza la hipótesis nula)
## 
## Prueba Dickey-Fuller para RESERVAS 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  series
## Dickey-Fuller = -2.7278, Lag order = 4, p-value = 0.279
## alternative hypothesis: stationary
## 
## RESERVAS es no estacionaria (no se rechaza la hipótesis nula)
## 
## Prueba Dickey-Fuller para POPULAR 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  series
## Dickey-Fuller = -2.5428, Lag order = 4, p-value = 0.3543
## alternative hypothesis: stationary
## 
## POPULAR es no estacionaria (no se rechaza la hipótesis nula)
## 
## Prueba Dickey-Fuller para JMMB_BDI 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  series
## Dickey-Fuller = -2.5589, Lag order = 4, p-value = 0.3478
## alternative hypothesis: stationary
## 
## JMMB_BDI es no estacionaria (no se rechaza la hipótesis nula)
#Vamos a crear un primer modelo para el Fondo Atlantico, para cada modelo los pasos serían: Separar los datos de entrenamiento, crear modelo, verificar los residuos (y si es necesario ajustarlo), pronosticar los datos y comparar los datos.

trainAtl <- window(ATLANTICO, start = c(2018, 1), end = c(2022, 12))
testAtl <- window(ATLANTICO, start = c(2023, 1))
plot.ts(trainAtl)

#Primero usaremos la función nntar en automatico, veremos que modelo sugiere la función y apartir de ahí veremos si lo modificamos

modeloAlt <- nnetar(trainAtl)
checkresiduals(modeloAlt)

## 
##  Ljung-Box test
## 
## data:  Residuals from NNAR(2,1,2)[12]
## Q* = 14.702, df = 12, p-value = 0.2581
## 
## Model df: 0.   Total lags used: 12
modAlt <- forecast(modeloAlt, h=12, level = 95)
autoplot(modAlt)

ProyeccionesAtlantico <- cbind(testAtl,modAlt$mean)
autoplot(ProyeccionesAtlantico)

rmse <- sqrt(mean((modAlt$mean - testAtl)^2))
mae <- mean(abs(modAlt$mean - testAtl))
mape <- mean(abs((modAlt$mean - testAtl) / testAtl)) * 100

cat("RMSE:", rmse, "\n")
## RMSE: 3.606821
cat("MAE:", mae, "\n")
## MAE: 3.308294
cat("MAPE:", mape, "%\n")
## MAPE: 30.93879 %
#En el modelo experimental, nos sugurió un (2,1,2), es decir dos retroceso estacionarios, uno estaciorario y 2 capas ocultas, para un periodo de 12, es decir mensual, vemos como resultado un MAPE de 34.85, RMSE 3.91 y MAE de 3.66. A continuación vamos a ajustar el modelo a que tenga 6 retroceso estacionarios, 2 no estacionarios y 5 capas oculta con 50 repeticiones.
modeloAlt <- nnetar(trainAtl,  p=5, P=2, size=5, repeats=50)
checkresiduals(modeloAlt)

## 
##  Ljung-Box test
## 
## data:  Residuals from NNAR(5,2,5)[12]
## Q* = 21.743, df = 12, p-value = 0.0405
## 
## Model df: 0.   Total lags used: 12
modAlt <- forecast(modeloAlt, h=12, level = 95)
autoplot(modAlt)

ProyeccionesAtlantico <- cbind(testAtl,modAlt$mean)
autoplot(ProyeccionesAtlantico)

rmse <- sqrt(mean((modAlt$mean - testAtl)^2))
mae <- mean(abs(modAlt$mean - testAtl))
mape <- mean(abs((modAlt$mean - testAtl) / testAtl)) * 100

cat("RMSE:", rmse, "\n")
## RMSE: 0.9467477
cat("MAE:", mae, "\n")
## MAE: 0.8325252
cat("MAPE:", mape, "%\n")
## MAPE: 8.073379 %
#En el modelo ajustado (6,2,5), vemos una mejora con un MAPE de 14.62, RMSE 1.97 y MAE de 14.62, para el caso de Afp Atlantico este es el mejor modelo neuronal, vamos a ir probando con los otros fondos

plot(modAlt)
lines(testAtl, col = "red")
legend("topright", legend = c("Datos de prueba","Predicciones"), col = c("red","blue"), lty = 3)

#MODELO CRECER

trainCre <- window(CRECER, start = c(2018, 1), end = c(2022, 12))
testCre <- window(CRECER, start = c(2023, 1))
plot.ts(trainCre)

modeloCre <- nnetar(trainCre, p=5, P=2, size=5, repeats=50)
checkresiduals(modeloCre)

## 
##  Ljung-Box test
## 
## data:  Residuals from NNAR(5,2,5)[12]
## Q* = 17.56, df = 12, p-value = 0.1297
## 
## Model df: 0.   Total lags used: 12
modCre <- forecast(modeloCre, h=12, level = 95)
autoplot(modCre)

ProyeccionesCrecer <- cbind(testCre,modCre$mean)
autoplot(ProyeccionesCrecer)

rmseCre <- sqrt(mean((modCre$mean - testCre)^2))
maeCre <- mean(abs(modCre$mean - testCre))
mapeCre <- mean(abs((modCre$mean - testCre) / testCre)) * 100

cat("RMSE:", rmseCre, "\n")
## RMSE: 1.29744
cat("MAE:", maeCre, "\n")
## MAE: 1.200054
cat("MAPE:", mapeCre, "%\n")
## MAPE: 12.13828 %
plot(modCre)
lines(testCre, col = "red")
legend("topright", legend = c("Datos de prueba","Predicciones"), col = c("red","blue"), lty = 1)

#MODELO JMMB_BDI

trainJmmb <- window(JMMB_BDI, start = c(2018, 1), end = c(2022, 12))
testJmmb <- window(JMMB_BDI, start = c(2023, 1))
plot.ts(trainJmmb)

modeloJmmb <- nnetar(trainJmmb,p=5, P=2, size=5, repeats=50)
checkresiduals(modeloJmmb)

## 
##  Ljung-Box test
## 
## data:  Residuals from NNAR(5,2,5)[12]
## Q* = 13.326, df = 12, p-value = 0.3458
## 
## Model df: 0.   Total lags used: 12
modJmmb <- forecast(modeloJmmb, h=12, level = 95)
plot(modJmmb)
lines(testJmmb, col = "red")
legend("topright", legend = c("Datos de prueba","Predicciones"), col = c("red","blue"), lty = 1)

rmseJmmb <- sqrt(mean((modJmmb$mean - testJmmb)^2))
maeJmmb <- mean(abs(modJmmb$mean - testJmmb))
mapeJmmb <- mean(abs((modJmmb$mean - testJmmb) / testJmmb)) * 100

cat("RMSE:", rmseJmmb, "\n")
## RMSE: 1.408177
cat("MAE:", maeJmmb, "\n")
## MAE: 1.202764
cat("MAPE:", mapeJmmb, "%\n")
## MAPE: 13.22916 %
#MODELO POPULAR
trainPop <- window(POPULAR, start = c(2018, 1), end = c(2022, 12))
testPop <- window(POPULAR, start = c(2023, 1))
plot.ts(trainPop)

modeloPop <- nnetar(trainPop, p=5, P=2,size=5, Repeats = 50)
checkresiduals(modeloPop)

## 
##  Ljung-Box test
## 
## data:  Residuals from NNAR(5,2,5)[12]
## Q* = 11.783, df = 12, p-value = 0.4632
## 
## Model df: 0.   Total lags used: 12
modPop <- forecast(modeloPop, h=12, level = 95)
plot(modPop)
lines(testPop, col = "red")
legend("topright", legend = c("Datos de prueba","Predicciones"), col = c("red","blue"), lty = 1)

rmsePop <- sqrt(mean((modPop$mean - testPop)^2))
maePop <- mean(abs(modPop$mean - testPop))
mapePop <- mean(abs((modPop$mean - testPop) / testPop)) * 100

cat("RMSE:", rmsePop, "\n")
## RMSE: 0.8283742
cat("MAE:", maePop, "\n")
## MAE: 0.7950945
cat("MAPE:", mapePop, "%\n")
## MAPE: 8.147243 %
#MODELO RESERVAS
trainRes <- window(RESERVAS, start = c(2018, 1), end = c(2022, 12))
testRes <- window(RESERVAS, start = c(2023, 1))
plot.ts(trainRes)

modeloRes <- nnetar(trainRes, p=5, P=2,size=5, Repeats = 50)
checkresiduals(modeloRes)

## 
##  Ljung-Box test
## 
## data:  Residuals from NNAR(5,2,5)[12]
## Q* = 10.387, df = 12, p-value = 0.5821
## 
## Model df: 0.   Total lags used: 12
modRes <- forecast(modeloRes, h=12, level = 95)
plot(modRes)
lines(testRes, col = "red")
legend("topright", legend = c("Datos de prueba","Predicciones"), col = c("red","blue"), lty = 1)

rmseRes <- sqrt(mean((modRes$mean - testRes)^2))
maeRes <- mean(abs(modRes$mean - testRes))
mapeRes <- mean(abs((modRes$mean - testRes) / testRes)) * 100

cat("RMSE:", rmseRes, "\n")
## RMSE: 0.6974309
cat("MAE:", maeRes, "\n")
## MAE: 0.5893661
cat("MAPE:", mapeRes, "%\n")
## MAPE: 6.282174 %
#MODELO ROMANA
trainRom <- window(ROMANA, start = c(2018, 1), end = c(2022, 12))
testRom <- window(ROMANA, start = c(2023, 1))
plot.ts(trainRom)

modeloRom <- nnetar(trainRom, p=5, P=2,size=5, Repeats = 50)
checkresiduals(modeloRom)

## 
##  Ljung-Box test
## 
## data:  Residuals from NNAR(5,2,5)[12]
## Q* = 31.263, df = 12, p-value = 0.001796
## 
## Model df: 0.   Total lags used: 12
modRom <- forecast(modeloRom, h=12, level = 95)
plot(modRom)
lines(testRom, col = "red")
legend("topright", legend = c("Datos de prueba","Predicciones"), col = c("red","blue"), lty = 1)

rmseRom <- sqrt(mean((modRom$mean - testRom)^2))
maeRom <- mean(abs(modRom$mean - testRom))
mapeRom <- mean(abs((modRom$mean - testRom) / testRom)) * 100

cat("RMSE:", rmseRom, "\n")
## RMSE: 1.624217
cat("MAE:", maeRom, "\n")
## MAE: 1.302793
cat("MAPE:", mapeRom, "%\n")
## MAPE: 12.49214 %
#MODELO SIEMBRA
trainSiem <- window(SIEMBRA, start = c(2018, 1), end = c(2022, 12))
testSiem <- window(SIEMBRA, start = c(2023, 1))
plot.ts(trainSiem)

modeloSiem <- nnetar(trainSiem, p=5, P=2,size=5, Repeats = 50)
checkresiduals(modeloSiem)

## 
##  Ljung-Box test
## 
## data:  Residuals from NNAR(5,2,5)[12]
## Q* = 11.575, df = 12, p-value = 0.4804
## 
## Model df: 0.   Total lags used: 12
modSiem <- forecast(modeloSiem, h=12, level = 95)
plot(modSiem)
lines(testSiem, col = "red")
legend("topright", legend = c("Datos de prueba","Predicciones"), col = c("red","blue"), lty = 1)

rmseSiem <- sqrt(mean((modSiem$mean - testSiem)^2))
maeSiem <- mean(abs(modSiem$mean - testSiem))
mapeSiem <- mean(abs((modSiem$mean - testSiem) / testSiem)) * 100

cat("RMSE:", rmseSiem, "\n")
## RMSE: 0.5435263
cat("MAE:", maeSiem, "\n")
## MAE: 0.4822435
cat("MAPE:", mapeSiem, "%\n")
## MAPE: 4.935217 %
#Una función para guardar los graficos y luego combinarlos

plot_function <- function(mod, test, titulo) {
  plot(mod, main = titulo, ylab = "Rentabilidad")
  lines(test, col = "red")
  legend("topright", legend = c("Datos de prueba", "Predicciones"), col = c("red", "blue"), lty = 3)
}

# Usar la función con diferentes modelos y datos de prueba

grafico_modAlt <- plot_function(modAlt, testAtl, "Modelo Atlántico")

grafico_modCre <- plot_function(modCre, testCre, "Modelo Crecer")

grafico_modJmmb <- plot_function(modJmmb, testJmmb, "Modelo Jmmb-BDI")

grafico_modPop <- plot_function(modPop, testPop, "Modelo Popular")

grafico_modRes <- plot_function(modRes, testRes, "Modelo Reservas")

grafico_modRom <- plot_function(modRom, testRom, "Modelo Romana")

grafico_modSiem <- plot_function(modSiem, testSiem, "Modelo Siembra")

library(cowplot)
## Warning: package 'cowplot' was built under R version 4.3.3
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:lubridate':
## 
##     stamp
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 4.3.3
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
library(forecast)
library(ggplot2)
library(cowplot)


# Convertir las predicciones en gráficos ggplot
grafico1 <- autoplot(modAlt, ylab = "Rentabilidad") + labs(title = "Modelo Atlántico")
grafico2 <- autoplot(modRom, ylab = "Rentabilidad") + labs(title = "Modelo Romana")
grafico3 <- autoplot(modCre, ylab = "Rentabilidad") + labs(title = "Modelo Crecer")
grafico4 <- autoplot(modRes, ylab = "Rentabilidad") + labs(title = "Modelo Reservas")
grafico5 <- autoplot(modJmmb, ylab = "Rentabilidad") + labs(title = "Modelo Jmmb-BDI")
grafico6 <- autoplot(modPop, ylab = "Rentabilidad") + labs(title = "Modelo Popular")
grafico7 <- autoplot(modSiem, ylab = "Rentabilidad") + labs(title = "Modelo Siembra")
# Combinar los gráficos en mosaicos separados
plot_grid(grafico1, grafico2,grafico3,grafico4,grafico5,grafico6, grafico7, nrow = 3)

#Combinando los graficos con sus datos de pruebas y datos pronosticados

grafico_prueba <- autoplot(testAtl) + labs(title = "Prueba")
datos_prueba <- ggplot_build(grafico_prueba)$data[[1]]
grafico1 <- grafico1 + geom_line(data = datos_prueba, aes(x = x, y = y), color ="red")

grafico_prueba2 <- autoplot(testRom) + labs(title = "Prueba")
datos_prueba2 <- ggplot_build(grafico_prueba2)$data[[1]]
grafico2 <- grafico2 + geom_line(data = datos_prueba2, aes(x = x, y = y), color ="red")

grafico_prueba3 <- autoplot(testCre) + labs(title = "Prueba")
datos_prueba3 <- ggplot_build(grafico_prueba3)$data[[1]]
grafico3 <- grafico3 + geom_line(data = datos_prueba3, aes(x = x, y = y), color ="red")

grafico_prueba4 <- autoplot(testRes) + labs(title = "Prueba")
datos_prueba4 <- ggplot_build(grafico_prueba4)$data[[1]]
grafico4 <- grafico4 + geom_line(data = datos_prueba4, aes(x = x, y = y), color ="red")

grafico_prueba5 <- autoplot(testJmmb) + labs(title = "Prueba")
datos_prueba5 <- ggplot_build(grafico_prueba5)$data[[1]]
grafico5 <- grafico5 + geom_line(data = datos_prueba5, aes(x = x, y = y), color ="red")

grafico_prueba6 <- autoplot(testPop) + labs(title = "Prueba")
datos_prueba6 <- ggplot_build(grafico_prueba6)$data[[1]]
grafico6 <- grafico6 + geom_line(data = datos_prueba6, aes(x = x, y = y), color ="red")

grafico_prueba7 <- autoplot(testSiem) + labs(title = "Prueba")
datos_prueba7 <- ggplot_build(grafico_prueba7)$data[[1]]
grafico7 <- grafico7 + geom_line(data = datos_prueba7, aes(x = x, y = y), color ="red")


# Combinar los gráficos en mosaicos separados con pronosticos vs pruebas

plot_grid(grafico1, grafico2,grafico3,grafico4,grafico5,grafico6, grafico7, nrow = 3)