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)
