# Instalar y cargar librerias necesarias para el proceso

library(fpp2)
library(readxl)
library(forecast)# Contiene el modelo ARIMA
library(tseries) #Para series de tiempo
library(TSA)     #Para series de tiempo
library(ggplot2) #Para hacer gráficos
library(dplyr)   #Para la manipulación de datos (filtrar, seleccionar, agregar, transformar)
library(stats)   #Se usa para diversas pruebas estadísticas (medias,varianza, arima,etc)
library(seasonal)#Para calcular la serie ajustada de estacionalidad
library(zoo)     #Para calcular la serie ajustada de estacionalidad
library(tidyverse)
library(lubridate)
library(tidyverse)
library(urca)

REDES NEURONALES RECURRENTES (Modelos Elman y Jordan)

Una red neuronal recurrente no tiene una estructura de capas definida, sino que permiten conexiones arbitrarias entre las neuronas, incluso pudiendo crear ciclos, con esto se consigue crear la temporalidad, permitiendo que la red tenga memoria.

Los RNN se denominan recurrentes porque realizan la misma tarea para cada elemento de una secuencia, y la salida depende de los cálculos anteriores.

MODELO ELMAN

En las redes de Elman, las entradas de estas neuronas, se toman desde las salidas de las neuronas de una de las capas ocultas, y sus salidas se conectan de nuevo en las entradas de esta misma capa, lo que proporciona una especie de memoria sobre el estado anterior de dicha capa.

library(readxl)
veh<- read_excel("D:/Veh.xlsx")
Y=ts(veh, start = c(2014,1),frequency = 12)
Y
##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 2014 20115 23744 24075 26128 26865 22974 27650 27465 29528 31386 25700 40393
## 2015 21241 22871 24671 21863 22525 22476 26595 23208 24855 22412 21004 29546
## 2016 16740 20278 19737 20291 19464 20589 17526 22409 21378 19168 22384 33431
## 2017 17215 17931 21022 17082 19898 19826 18607 20901 19297 19689 21103 25386
## 2018 16399 18293 19523 20166 20522 18270 19315 21410 21370 21496 26245 33043
## 2019 15941 18440 20170 19763 22269 19455 22874 23284 22659 23858 23946 30661
## 2020 18427 20547 12290   217  8933 11981 14435 13209 18408 20858 22330 26854
## 2021 14327 19648 22914 19021 14684 20417 23117 21052 22824 23340 23616 25312
## 2022 17391 19199 20837 20622 22411 23306 23233 24386 23871 22577 22625 21880
## 2023 13852 15761
#Damos formato de serie de tiempo
Y1=ts(veh, start = c(2014,1), end=c(2022,02),frequency = 12)
length(Y1)
## [1] 98
#Para aplicar redes neuronales debemos normalizar datos
#Para ello hemos asociado a nuestro dataset de base una variable “Z” y a partir de esta hemos realizar la normalización a través de la variable “S”.
Z <- as.ts(Y1,F)
S <- (Z-min(Z))/(max(Z)-min(Z))  
plot(S)

#A continuación comprobamos el numero de filas totales que contiene nuestro dataset.

tamano_total <- length(S)
tamano_total
## [1] 98
#Ahora dividiremos los conjuntos de entrenamiento en un 75% y prueba en un 25% respectivamente.
tamano_train <- round(tamano_total*0.75, digits = 0)
train <- 0:(tamano_train-1)
train
##  [1]  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
## [26] 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
## [51] 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73
test <- (tamano_train):tamano_total
test
##  [1] 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98
#Ahora crearemos un dataframe con n columnas, cada una de las cuales adelantara 
#un valor de la serie en el futuro, a través de una variable tipo zoo, equivalente al #periodo de retardo de la serie.

library(zoo)
library(quantmod)
y <- as.zoo(S)
x1 <- Lag(y, k = 1)
x2 <- Lag(y, k = 2)
x3 <- Lag(y, k = 3)
x4 <- Lag(y, k = 4)
x5 <- Lag(y, k = 5)
x6 <- Lag(y, k = 6)
x7 <- Lag(y, k = 7)
x8 <- Lag(y, k = 8)
x9 <- Lag(y, k = 9)
x10 <- Lag(y, k = 10)
x11 <- Lag(y, k = 11)
x12 <- Lag(y, k = 12)
slogN <- cbind(y,x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12)
DT::datatable(slogN)
#A continuacion eliminaremos los valores NA producidos al desplazar la serie:
slogN1 <- slogN[-(1:12),]
DT::datatable(slogN1)
#Luego definimos los valores de entrada y salida de la red neuronal:
inputs <- slogN1[,2:13]
outputs <- slogN1[,1]
#Ahora crearemos la red de Elman, probando diferentes tipos de combinaciones de neuronas en las capas ocultas e iteraciones máximas, ademas del ritmo de aprendizaje, para ajustar lo mejor posible la curva de predicción a la del modelo de la serie. De esta forma hemos llegado a estos valores a la hora de crear nuestra red. Asi mismo ponemos una semilla para que el resultado sea reproducible.

library(RSNNS)
## Loading required package: Rcpp
set.seed(42)
fit<-elman(inputs[train],outputs[train],size=5,learnFuncParams=c(0.1),
                  maxit=10000)
#En la gráfica siguiente vemos como evoluciona el error de la red con el numero de iteraciones para los parámetros expuestos.

plotIterativeError(fit, main = "Iterative Error for 5 Neuron")

En la gráfica anterior se observa una convergencia rápida hacia cero, se espera entonces un ajuste rápido del modelo.

#Ahora realizamos la predicción con el resto de los términos de la serie que son los datos #seleccionados para test, pasamos pues una vez entrenada a probarla y a representarla #graficamente para ver el ajuste del modelo.

y <- as.vector(outputs[-test])
plot(y,type="l")
pred <- predict(fit, inputs[-test])
lines(pred,col = "red")

El ajuste predice bastante bien con los parametros elegidos, pues la curva del modelo de la serie y la de la prediccion parecen bastante ajustadas.

Esta representacion grafica se puede utilizar para ir ajustando la prediccion y el modelo a medida que vamos probando diferentes parametros de la red de Elman, de forma que la curva del modelo y de la prediccion queden lo mas ajustados posibles.

#Ahora gracias al efecto memoria vamos a adelantarle a la serie al menos en un valor con una precision muy buena. Para ello volveremos a introducir los datos de entrenamiento.

predictions <- predict(fit,inputs[-train])
predictions
##               [,1]
## feb 2021 0.3329012
## mar 2021 0.2609815
## abr 2021 0.0670051
## may 2021 0.4164404
## jun 2021 0.3693653
## jul 2021 0.4730132
## ago 2021 0.4960409
## sep 2021 0.5433677
## oct 2021 0.5297236
## nov 2021 0.6084241
## dic 2021 0.6646473
## ene 2022 0.3588158
## feb 2022 0.5120363
#posteriori desnormalizaremos los datos:

mod1 <- predictions*(max(Z)-min(Z))+min(Z)
mod1
##               [,1]
## feb 2021 13591.640
## mar 2021 10702.193
## abr 2021  2908.997
## may 2021 16947.908
## jun 2021 15056.622
## jul 2021 19220.779
## ago 2021 20145.938
## sep 2021 22047.342
## oct 2021 21499.177
## nov 2021 24661.048
## dic 2021 26919.869
## ene 2022 14632.786
## feb 2022 20788.571

Aquí vemos la gráfica con los valores pronosticados con la linea roja. -Los valores que adelantamos en el tiempo corresponden a mod1, de los cuales adelantaremos 12 meses a futuro para nuestro estudio.

#Ahora veamos la representación de los valores predecidos para el siguiente periodo.
x <- 1:(tamano_total+length(mod1))
y <- c(as.vector(Z),mod1)
plot(x[1:tamano_total], y[1:tamano_total],col = "blue", type="l")
lines( x[(tamano_total):length(x)], y[(tamano_total):length(x)], col="red")

length(y)
## [1] 111

MODELO JORDAN

En las redes Jordan, la diferencia esta en que la entrada de las neuronas de la capa de contexto se toma desde la salida de la red.

Realizamos las mismas operaciones que con la red Elman, sustituyendo el modelo, obtenemos el resultado para la red Jordan.

set.seed(42)
fit <-jordan(inputs[train],outputs[train],size=5,learnFuncParams=c(0.1),
             maxit=35000)

plotIterativeError(fit, main = "Iterative Error for 5 Neuron")

y <- as.vector(outputs[-test])
plot(y,type="l")
pred <- predict(fit, inputs[-test])
lines(pred,col = "red")

predictions <- predict(fit,inputs[-train])
mod2 <- predictions*(max(Z)-min(Z))+min(Z)
mod2
##               [,1]
## feb 2021  3884.444
## mar 2021 33576.846
## abr 2021 11833.495
## may 2021 14823.377
## jun 2021 13574.293
## jul 2021 20448.654
## ago 2021 18429.020
## sep 2021 19670.796
## oct 2021 20747.208
## nov 2021 27825.488
## dic 2021 25797.013
## ene 2022 13655.142
## feb 2022 17335.238
x <- 1:(tamano_total+length(mod2))
y <- c(as.vector(Z),mod2)
plot(x[1:tamano_total], y[1:tamano_total],col = "blue", type="l")
lines( x[(tamano_total):length(x)], y[(tamano_total):length(x)], col="red")

Estimación del error Comparativo de los modelos con los valores actuales observados.

data=veh[99:110,]
data_real <- ts(data, start = c(2022,3), end=c(2023,02), frequency = 12)

data_real
##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 2022             20837 20622 22411 23306 23233 24386 23871 22577 22625 21880
## 2023 13852 15761
#Para Elman y Jordan al no ser modelos con forecast, lo convertimos en series de tiempo #para que lo acepte el comando accuracy.De tal forma que:

mod1[1:12]
##  [1] 13591.640 10702.193  2908.997 16947.908 15056.622 19220.779 20145.938
##  [8] 22047.342 21499.177 24661.048 26919.869 14632.786
m1 <- mod1[1:12]
mod1c <- ts(m1, frequency=12,start=c(2022,3))
mod1c
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2022                     13591.640 10702.193  2908.997 16947.908 15056.622
## 2023 26919.869 14632.786                                                  
##            Aug       Sep       Oct       Nov       Dec
## 2022 19220.779 20145.938 22047.342 21499.177 24661.048
## 2023
m2 <- mod2[1:12]
mod2c <- ts(m2, frequency=12,start=c(2022,3))
mod2c
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2022                      3884.444 33576.846 11833.495 14823.377 13574.293
## 2023 25797.013 13655.142                                                  
##            Aug       Sep       Oct       Nov       Dec
## 2022 20448.654 18429.020 19670.796 20747.208 27825.488
## 2023

Elman (RRN)

accuracy(mod1c,data_real)
##                ME     RMSE      MAE      MPE     MAPE      ACF1 Theil's U
## Test set 3918.892 8467.508 6560.378 14.71541 32.55703 0.5280143   3.23668

Jordan

accuracy(mod2c,data_real)
##                ME     RMSE     MAE      MPE    MAPE       ACF1 Theil's U
## Test set 2591.269 8999.551 7732.16 8.647962 38.0191 -0.1205701   2.92748

De acuerdo a las métricas de evaluación, el mejor ajuste en predicción entre los dos modelos, lo arroja el modelo de red recurrente “ELMAN”, MAE: 6.560 y RMSE: 8467.