0. EVOLUCION IPC ESPAÑA.

-Hemos escogido para este estudio un extracto de datos referente a la “Evolucion mensual del Indice de Precios al Consumo (IPC) en España”, en el periodo entre Enero del 2012 y Enero del 2019.

Fuente: https://www.idescat.cat/indicadors/?id=conj&n=10261&lang=es&col=1

-De este mismo extracto vamos a pronosticar el indice general de precios al consumo de los siguientes 8 meses despues, para ello utilizaremos como valores reales para comprobar y comparar con los pronosticos,los datos de esta misma fuente desde Febrero del 2019 a Septiembre del 2019.

-Una vez hecha las predicciones y comparados todos los modelos con los resultados reales, procederemos a estimar el error de cada modelo, para posteriormente elegir el mejor o los mejores modelos de pronostico de todos los que hemos representado.

1. MODELO ARIMA.

-Un modelo autorregresivo integrado de promedio movil o ARIMA es un modelo estadistico que utiliza variaciones y regresiones de datos estadisticos con el fin de encontrar patrones para una prediccion hacia el futuro. Se trata de un modelo dinamico de series temporales, es decir, las estimaciones futuras vienen explicadas por los datos del pasado y no por variables independientes.

CUESTIONES IMPORTANTES A TENER EN CUENTA

-Para definir un modelo ARIMA, debemos tener en cuenta que los datos deben ser estacionarios,es decir cuando su distribucion y sus parametros no varian con el tiempo.

-Otra cuestion para tener en cuenta es que los datos deben ser univariantes, ARIMA trabaja en una sola variable.

La regresion automatica tiene que ver con los valores pasados, es decir necesitamos un historico de los datos para poder realizar el modelo y asi predecir datos futuros.

PROCESO DEL MODELO

-Comenzamos pues cargando las librerias previamente instaladas, que nos hacen falta para ejecutar el modelo tal como “tidyverse” para realizar transformaciones en general de nuestro conjunto de datos, ademas de otras funciones asociadas , “readxl” para leer los datos de un archivo excel, “lubridate” para hacer transformaciones en un campo fecha, y forecast que es una libreria empleada para pronosticos de series de tiempo, asi como tseries

library(tidyverse)
library(readxl)
library(forecast)
library(lubridate)
library(tseries)

-A continuacion vamos a cargar nuestros conjunto de datos del archivo excel que los contiene:

p_ipc <- read_excel("Evolucion_IPC_es.xlsx")
p_ipc

## # A tibble: 85 x 14
##    Periodo             `Índice general` `Alimentos y bebida~ `Bebidas alcohólic~
##    <dttm>                         <dbl>                <dbl>               <dbl>
##  1 2012-01-01 00:00:00             97.4                 94.3                88.1
##  2 2012-02-01 00:00:00             97.5                 94.4                88.4
##  3 2012-03-01 00:00:00             98.1                 94.5                88.4
##  4 2012-04-01 00:00:00             99.5                 94.7                91  
##  5 2012-05-01 00:00:00             99.3                 94.5                91.9
##  6 2012-06-01 00:00:00             99.2                 94.9                91.9
##  7 2012-07-01 00:00:00             98.9                 94.8                92.2
##  8 2012-08-01 00:00:00             99.5                 95.2                92.3
##  9 2012-09-01 00:00:00            100.                  95.6                92.6
## 10 2012-10-01 00:00:00            101.                  95.9                92.6
## # ... with 75 more rows, and 10 more variables: Vestido y calzado <dbl>,
## #   Vivienda <dbl>, Menaje del hogar <dbl>, Medicina <dbl>, Transportes <dbl>,
## #   Comunicaciones <dbl>, Ocio y cultura <dbl>, Enseñanza <dbl>,
## #   Hotelescafés y restaurantes <dbl>, Otros bienes y servicios <dbl>

Vemos que nuestro dataset esta compuesto por diferentes columnas donde se encuentran el “Periodo” o Fecha(year/month/day), otra el “Índice general” que es la variable a predecir y las demas columnas constituidas por los distintos indices de precios establecidos segun sectores o productos.
Ahora vamos a hacer las transformaciones dentro del dataset en el campo “Periodo” correspondientes para separar month y year y luego añadir la variable objetivo de pronostico “Índice general”, para configurar nuestro conjunto de datos del cual vamos a partir para iniciar nuestro pronostico, todo ello guardado en la variable p_ipc_es, donde mostramos los resultados de la transformacion.

p_ipc %>%
        mutate_at(vars(Periodo), list(year=year, month=month, day=day))%>%
        select("month", "year","Índice general")->p_ipc_es

p_ipc_es

## # A tibble: 85 x 3
##    month  year `Índice general`
##    <dbl> <dbl>            <dbl>
##  1     1  2012             97.4
##  2     2  2012             97.5
##  3     3  2012             98.1
##  4     4  2012             99.5
##  5     5  2012             99.3
##  6     6  2012             99.2
##  7     7  2012             98.9
##  8     8  2012             99.5
##  9     9  2012            100. 
## 10    10  2012            101. 
## # ... with 75 more rows

Ahora vamos a transformar nuestro conjunto de datos en una serie temporal que asignaremos a una variable Y, indicando el periodo y la frecuencia = 12 (meses), a la que vamos a someter a estudio. La variable objetivo a pronosticar sera el “Índice general”

-El resultado nos dara la transformacion de nuestros datos originales a formato “ts”, que nos indica que es una serie temporal.

Y <- ts(p_ipc_es[,3], start = c(2012,1), end = c(2019,1), frequency = 12)

Y

##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 2012  97.4  97.5  98.1  99.5  99.3  99.2  98.9  99.5 100.5 101.3 101.2 101.3
## 2013 100.0 100.1 100.5 100.9 101.0 101.2 100.7 101.0 100.8 101.2 101.4 101.5
## 2014 100.2 100.1 100.3 101.2 101.3 101.3 100.3 100.5 100.6 101.1 101.0 100.5
## 2015  98.8  99.0  99.7 100.6 101.1 101.3 100.4 100.1  99.8 100.4 100.8 100.5
## 2016  98.6  98.2  98.8  99.5 100.1 100.5  99.8  99.9  99.9 101.1 101.4 102.0
## 2017 101.5 101.1 101.1 102.1 102.0 102.1 101.4 101.6 101.7 102.7 103.1 103.2
## 2018 102.1 102.2 102.3 103.2 104.1 104.4 103.6 103.8 104.0 105.0 104.9 104.4
## 2019 103.1

class(Y)

## [1] "ts"

Ahora vamos a dar una visual a los datos que nos ofrece la serie, graficamente.

autoplot(Y) +
  ggtitle("IPC Consumo en España") +
  ylab("Índice General")

A priori damos una visual a los datos y apreciamos que la serie parece tener ciertos patrones comunes.
Vamos a ver ahora un grafico de la descomposicion de la serie la cual se descompone en 3 graficos: seasonal, trend y el grafico donde se representan los residuos o remainder.

descomp = decompose(Y)
autoplot(descomp)

En “data” podemos apreciar la grafica de los datos, en “seasonal” que es la estacionalidad, apreciamos el patron de los datos que se repite cada año, que es el patron que se mantiene constante.
Por otro lado “trend” representa la tendencia, que es como se puede ver tiene primero una pequeña subida luego baja y luego a partir de ese punto toma una tendencia ascendente, y “remainder”, representa la configuracion de los residuos que va dejando el modelo, cuando disminuyen y cuando aumentan. Si sumamos las graficas de seasonal, trend y remainder nos da la grafica total data de la serie temporal

-Vamos a comprobar ahora para asegurarnos si nuestra serie temporal es estacionaria o no a traves del test de Dickey-Fuller, si cumple en su ejecucion que el p-valor>0.05 la serie cumple la hipotesis nula en la cual la serie no es estacionaria, si por el contrario el p-valor<0.05 no se cumple la hipotesis nula con lo cual la serie es estacionaria, Asi pues lo comprobamos:

adf.test(Y)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  Y
## Dickey-Fuller = -1.401, Lag order = 4, p-value = 0.8218
## alternative hypothesis: stationary

Y vemos que el p-valor nos sale 0.8218>0.05 por tanto se cumple la hipotesis nula y la serie no es estacionaria.Asi que tendremos que buscar cambios para hacer que la serie se convierta en estacionaria y asi aplicar el modelo ARIMA.
Ahora vamos a utilizar 2 parametros “acf” y “pacf” que miden el nivel de residuos, de forma que los podemos ver graficamente. En la grafica aparece el nivel de los residuos y los rezagos(Lags), para identificar la serie de tiempo que sea autoregresiva los valores no deberian sobrepasar el nivel de significancia establecido representado por las lineas de rayas horizontales entrecortadas, como podemos ver los niveles de significancia son del 2% aproximadamente. En ambos casos hacen un pico de rezago que decrece hacia 0, donde observamos que en el caso de acf los residuos superan ampliamente los niveles de significancia, y en el de pacf solo se observa en alguno de sus niveles arriba y abajo por encima y por abajo del 0 que sobrepasan los niveles de significancia, ese 2%.

acf(Y)

pacf(Y)

Vamos a ver con el comando ndiffs(), que diferencias hay para poder arreglar nuestra serie temporal y poder convertirla a estacionaria y asi poder aplicar el modelo ARIMA.

ndiffs(Y)

## [1] 1

El comando ndiffs nos indica que hay una diferencia, que tenemos que arreglar, vamos a determinar la grafica que nos da cuando aplicamos la diferencia en el tiempo a la serie.

DY <- diff(Y, lag = 1)

autoplot(DY) +
  ggtitle("Diferencias IPC España") +
  ylab("Índice General")

La grafica nos muestra la diferencia que hemos aplicado a la serie, en la que se encuentra con una media en torno a 0 y una varianza entre 1 y -1 mas o menos de la mayoria de los datos, salvando algunos picos, por tanto, constatamos su estacionariedad.

-Si aplicamos ahora el test de Dickey-Fuller despues de la diferenciacion nos confirma y nos da un resultado < 0.05 lo que confirma que la serie ahora si es estacionaria.

adf.test(DY)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  DY
## Dickey-Fuller = -5.1291, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

-Pasamos a continuacion a elaborar el modelo de pronostico ARIMA, con la funcion auto.arima. la funcion auto.arima esta diseñada para una estimacion rapida del modelo ARIMA probando varias series de tiempo.

-La funcion auto.arima, realiza la funcion ARIMA de forma automatica, va a ir iterando, va a ir buscando el modelo que mas se ajusta a los datos, y nos da respuesta del mejor modelo, como puede verse al ejecutar el print del mismo. El resultado despues de sucesivas iteraciones, se queda con el modelo que tuvo el menor valor de AIC.

model <- auto.arima(DY, trace = T, stepwise = F, approximation =F, max.d = 0)

print(model)

## Series: DY 
## ARIMA(1,0,1)(0,1,1)[12] 
## 
## Coefficients:
##           ar1     ma1     sma1
##       -0.2969  0.7494  -0.8815
## s.e.   0.1888  0.1333   0.3868
## 
## sigma^2 = 0.09695:  log likelihood = -24.73
## AIC=57.47   AICc=58.06   BIC=66.57

-Hacemos ahora pues del modelo un chequeo de los resultados:

checkresiduals(model)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,1)(0,1,1)[12]
## Q* = 10.088, df = 14, p-value = 0.7557
## 
## Model df: 3.   Total lags used: 17

-Vemos en los residuos que los valores perdidos en general estan dentro del limite de significancia (grafica ACF).

-Tambien vemos que la distribucion de los residuos (residuals), tiende a centrarse en 0 y van teniendo una distribucion normal.

-Vamos a hacer ahora una estimacion de este modelo, haciendo un forecast en el que vamos a predecir el valor de las ventas para los proximos 8 meses con un nivel de significancia de un 95%.

mod1 <- forecast(model, 8, level = 95)

-El resultado de la prediccion que nos da, no son datos reales, sino que corresponden a la serie diferenciada, tendremos que transformar estos datos a datos reales, de manera que hallamos el resultado real sumando la suma acumulada de los valores predichos mas el ultimo valor del indice de la serie original, de tal manera que:

last_Y = Y[85]
pred_mod1 = mod1$mean
pron_real = cumsum(pred_mod1) + last_Y

-Ahora vamos a convertir el resultado que nos da a un objeto ts que empiece en el indice despues del ultimo de los datos de Y, de tal manera que:

pron_real_ts = ts(pron_real, start = c(2019,2), frequency=12)

-Con esto tendremos el pronostico para la serie predicha en formato “ts”.

-Si agrupamos la serie original mas la serie predicha nos dara el plot de la serie total resultante, que acontinuacion podremos graficar, de tal manera que:

ts.plot(Y, pron_real_ts, lty = c(1,3), col = c(5,2))

-Ahora vamos a aplicar la libreria dygrahs para un grafico mas vistoso e interactivo.

library(dygraphs)

graf_inter = cbind(Y,pron_real_ts)

dygraph(graf_inter)

2. MODELO SUAVIZADO EXPONENCIAL(ETS).

-El algoritmo ETS es especialmente util para conjuntos de datos con estacionalidad y otras suposiciones previas sobre los datos. ETS calcula un promedio ponderado sobre todas las observaciones en el conjunto de datos de las series temporales de entrada como su prediccion. Las ponderaciones disminuyen exponencialmente con el tiempo, en lugar de las ponderaciones constantes en los metodos de promedio movil simple. Las ponderaciones dependen de un parametro constante, conocido como parametro de suavizamiento.

-Ahora vamos a desarrollar con este dataset un modelo de suavizado exponencial usando la funcion ets con los parametros por defecto. Luego vamos a realizar la prediccion de los proximos 8 meses y representaremos graficamente nuestra prediccion, asi mismo visualizaremos los valores estimados y el grado de certeza del 95% en los que se muestra los intervalos inferior y superior de los mismos.

-En primer lugar vamos a aplicar el modelo y ver visualmente la representacion de los residuos del mismo y vemos que se encuentran practicamente todos dentro de los limites de significancia del 2%, a excepcion de uno de los niveles (grafica ACF), y en general tienden a centrarse en 0, y van teniendo una distribucion normal(grafica residuals).

fit_ets_default <- ets(Y)
checkresiduals(fit_ets_default)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(M,N,A)
## Q* = 17.17, df = 3, p-value = 0.0006521
## 
## Model df: 14.   Total lags used: 17

-En segundo lugar vamos a calcular el forecast de nuestro modelo en los proximos 8 meses con un nivel de significancia del 95% y vamos a representarlo graficamente.

mod2 <- forecast(fit_ets_default, 8, level = 95)
plot(mod2)

-Ahora vamos a ver representado el ajuste entre los datos de la serie y el pronostico del modelo en la siguiente representacion grafica, para ello utilizamos la funcion fitted() que obtiene un ajuste con la data historica.

autoplot(mod2)+autolayer(fitted(mod2), series="Ajuste")

-Visualizamos en la grafica que la serie del modelo en negrita se ajusta bastante bien a la prediccion linea roja, y su correspondiente residual que se encuentra entre ambas lineas.

-Y visualizamos el pronostico, con sus margnenes inferior y superior del 95% de valores de grado de certeza (area de color lila)

-Tambien imprimimos un resumen completo de los resultados.

print(summary(mod2))

## 
## Forecast method: ETS(M,N,A)
## 
## Model Information:
## ETS(M,N,A) 
## 
## Call:
##  ets(y = Y) 
## 
##   Smoothing parameters:
##     alpha = 0.9985 
##     gamma = 0.0015 
## 
##   Initial states:
##     l = 98.0592 
##     s = 0.6565 0.8035 0.6712 -0.072 -0.1801 -0.379
##            0.4183 0.3851 0.1607 -0.6357 -0.997 -0.8314
## 
##   sigma:  0.0034
## 
##      AIC     AICc      BIC 
## 208.1662 215.1228 244.8060 
## 
## Error measures:
##                      ME      RMSE       MAE       MPE      MAPE      MASE
## Training set 0.06916369 0.3084548 0.2472397 0.0681582 0.2449994 0.2247633
##                   ACF1
## Training set 0.2998942
## 
## Forecasts:
##          Point Forecast    Lo 95    Hi 95
## Feb 2019       102.9333 102.2567 103.6099
## Mar 2019       103.2936 102.3357 104.2514
## Apr 2019       104.0908 102.9143 105.2673
## May 2019       104.3147 102.9535 105.6760
## Jun 2019       104.3488 102.8249 105.8726
## Jul 2019       103.5512 101.8827 105.2197
## Aug 2019       103.7489 101.9468 105.5511
## Sep 2019       103.8573 101.9306 105.7841

-Y podemos representar solo los resultados a traves de un dataframe, asignando la ejecucion del mismo a la variable pronostico.

pronostico <- as.data.frame(mod2)
pronostico

##          Point Forecast    Lo 95    Hi 95
## Feb 2019       102.9333 102.2567 103.6099
## Mar 2019       103.2936 102.3357 104.2514
## Apr 2019       104.0908 102.9143 105.2673
## May 2019       104.3147 102.9535 105.6760
## Jun 2019       104.3488 102.8249 105.8726
## Jul 2019       103.5512 101.8827 105.2197
## Aug 2019       103.7489 101.9468 105.5511
## Sep 2019       103.8573 101.9306 105.7841

3. SUAVIZADO EXPONENCIAL(HOLT-WINTERS)

-El metodo se basa en un algoritmo iterativo que a cada tiempo realiza un pronostico sobre el comportamiento de la serie en base a promedios debidamente ponderados de los datos obtenidos anteriormente.

-A este particular hay que reseñar los 2 diferentes tipos de estacionalidad que se pueden dar en las graficas, que son estacionalidad aditiva o estacionalidad multiplicativa.

-Elijiremos el modelo multiplicativo cuando la magnitud del patron estacional en los datos depende de la magnitud de los datos. En otras palabras, la magnitud del patron estacional aumenta a medida que los valores de los datos se incrementan y disminuye a medida que los valores de los datos decrecen.

-Elijiremos el modelo aditivo cuando la magnitud del patron estacional en los datos no dependa de la magnitud de los datos. En otras palabras, la magnitud del patron estacional no cambia cuando la serie sube o baja.

-vemos ahora el plot de nuestros datos:

plot(Y)

-Como la serie de partida es no estacional como pudimos ver anteriormente en el modelo arima, vamos aplicar a holtwinters, seasonal = “multiplicative”

-Asi que haremos una grafica donde se vea la serie del modelo (linea negra) y el ajuste con la prediccion(linea roja), asi mismo aparecera la descomposicion de la grafica para evaluar los datos.

m1 = HoltWinters(Y, seasonal = "multiplicative")

## Warning in HoltWinters(Y, seasonal = "multiplicative"): optimization
## difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH

plot(m1)

plot(fitted(m1))

-A posterori realizaremos las predicciones a 8 meses vista y las graficaremos:

mod3=predict(m1, 8, prediction.interval = TRUE)
mod3

##               fit      upr      lwr
## Feb 2019 103.0017 103.8162 102.1873
## Mar 2019 103.4045 104.4936 102.3155
## Apr 2019 104.1534 105.4673 102.8395
## May 2019 104.3332 105.8377 102.8287
## Jun 2019 104.4485 106.1237 102.7734
## Jul 2019 103.5586 105.3766 101.7406
## Aug 2019 103.7922 105.7594 101.8249
## Sep 2019 104.3191 106.4322 102.2059

plot(m1, mod3)

Asi mismo parece que la curva se ajusta bastante bien al modelo, y se muestra la prediccion a 8 meses con los margenes del intervalo (mod3).

3. REDES NEURONALES DE RETROALIMENTACION (nntear)

Son redes con una sola capa oculta y entradas retrasadas para pronosticar series de tiempo univariadas.
Vamos a elaborar entonces el modelo de red neuronal con nuestros datos:

set.seed(42)
neural_network <- nnetar(Y)

-Una vez hecho esto vamos a comprobar los residuos que presenta nuestro modelo.

checkresiduals(neural_network)

## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.

-En este caso hay niveles de residuos que sobrepasan la linea azul del 2% de significancia (Grafica ACF). y en general tienden a centrarse en 0, y van teniendo una distribucion normal(grafica residuals).

Una vez hecho esto realizamos la prediccion a 8 meses de los datos, con una significancia del 95% de los mismos.

mod4 <- forecast(neural_network, h=8, level = 95)

-Luego realizamos un grafico de nuestro pronostico.

autoplot(mod4)

-Ahora vamos a ver representado el ajuste entre lo datos de la serie y el pronostico del modelo en la siguiente representacion grafica, para ello utilizamos la funcion fitted() que obtiene un ajuste con la data historica.

autoplot(mod4)+autolayer(fitted(mod4), series="Ajuste")

En esta grafica se puede ver el residual que seria la diferencia entre la linea del modelo en negrita y la linea roja su prediccion.
Ahora vamos a determinar los resultados de este modelo agrupandolos en un dataframe y, asignandolo a la variable Pronostico para su posterior ejecucion.

pronostico <- as.data.frame(mod4)
pronostico

##           Feb      Mar      Apr      May      Jun      Jul      Aug      Sep
## 2019 103.5525 103.8891 104.2772 104.5930 104.7720 104.7571 104.7808 104.8156

4. 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 calculos anteriores.

4.1 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.

-Para empezar a desarrollar este modelo junto con el de Jordan, vamos a necesitar un par de librerias mas:

require(RSNNS)
require(quantmod)

-Ahora vamos a desarrollar el proceso para pronosticar a 8 meses de nuestros datos con esta red.

-En primer lugar para actuar en dicho proceso con redes neuronales tenemos que normalizar nuestros datos para que tomen valores entre 0 y 1. Para ello hemos asociado a nuestro dataset de base una variable “Z” y a raiz de esta hemos realizado la normalizacion a traves de la variable “S” que a posteriori representamos para su visualizacion.

Z <- as.ts(Y,F)
S <- (Z-min(Z))/(max(Z)-min(Z))  
plot(S)

-A continuacion comprobamos el numero de filas totales que contiene nuestro dataset y dividiremos los conjuntos de entrenamiento en un 75% y prueba en un 25% respectivamente.

tamano_total <- length(S)
tamano_total

## [1] 85

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

test <- (tamano_train):tamano_total
test

##  [1] 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85

-Ahora crearemos un dataframe con n columnas, cada una de las cuales adelantara un valor de la serie en el futuro, a traves de una variable tipo zoo, equivalente al periodo de retardo de la serie.

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)

-A continuacion eliminaremos los valores NA producidos al desplazar la serie:

slogN <- slogN[-(1:12),]

-Luego definimos los valores de entrada y salida de la red neuronal:

inputs <- slogN[,2:13]
outputs <- slogN[,1]

-Ahora crearemos la red de elman, probando diferentes tipos de combinaciones de neuronas en las capas ocultas e iteraciones maximas, ademas del ritmo de aprendizaje aunque este ultimo apenas lo hemos tocado, para ajustar lo mejor posible la curva de prediccion 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.

set.seed(42)
fit <- fit<-elman(inputs[train],outputs[train],size=c(7,3),learnFuncParams=c(0.1),maxit=64000)

-En la grafica siguiente vemos como evoluciona el error de la red con el numero de iteraciones para los parmetros expuestos.

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

-y vemos que el error converge a 0 muy rapidamente.

-Ahora pasamos a hacer la prediccion con el resto de los terminos 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")

-y parece que hemos realizado un ajuste que 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])

-y a posteriori desnormalizaremos los datos:

mod5 <- predictions*(max(Z)-min(Z))+min(Z)
mod5

##               [,1]
## abr. 2018 103.7537
## may. 2018 103.4365
## jun. 2018 104.1398
## jul. 2018 102.6478
## ago. 2018 103.2958
## sep. 2018 103.5337
## oct. 2018 105.0432
## nov. 2018 105.0289
## dic. 2018 105.0369
## ene. 2019 104.2651

-Ahora veremos la representacion de los valores predecidos para el siguiente periodo.

x <- 1:(tamano_total+length(mod5))
y <- c(as.vector(Z),mod5)
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")

-Aqui vemos la grafica con los valores predecidos con la linea roja.

-Los valores que adelantamos en el tiempo corresponden a mod5, de los cuales adelantaremos 8 meses a futuro para nuestro estudio.

4.2 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=6,learnFuncParams=c(0.1),maxit=78000)

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

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

Mas o menos la curva del modelo y de la prediccion se encuentran ajustadas.

predictions <- predict(fit,inputs[-train])

mod6 <- predictions*(max(Z)-min(Z))+min(Z)
mod6

##               [,1]
## abr. 2018 104.0573
## may. 2018 106.0169
## jun. 2018 106.6275
## jul. 2018 103.9661
## ago. 2018 103.2639
## sep. 2018 102.8324
## oct. 2018 104.5126
## nov. 2018 109.1480
## dic. 2018 108.3287
## ene. 2019 106.7435

x <- 1:(tamano_total+length(mod6))
y <- c(as.vector(Z),mod6)
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")

-Aqui vemos la grafica con los valores predecidos con la linea roja.

-Los valores que adelantamos en el tiempo corresponden a mod6, de los cuales adelantaremos 8 meses a futuro para nuestro estudio.

4.3 ESTIMACION DEL ERROR COMPARATIVO DE LOS MODELOS CON LOS VALORES ACTUALES OBSERVADOS.

-Estos son los valores actuales observados (datos reales pasado el tiempo de haber desarrollado el pronostico), comparado con los valores predecidos por cada uno de los modelos sometidos a estudio:

                Periodo Observado    Arima     ETS     Holt-W    nntear   Elman(RNN)  Jordan(RNN) 

                feb.-19  103.3       103.0    102.9    103.0       103.4       103.7          104.0
                mar.-19  103.7       103.4    103.2    103.4       103.7       103.4          106.0
                abr.-19  104.7       104.3    104.0    104.1       104.0       104.1          106.6
                may.-19  104.9       104.6    104.3    104.3       104.2       102.6          103.9
                jun.-19  104.8       104.7    104.3    104.4       104.4       103.2          103.2
                jul.-19  104.2       104.0    103.5    103.5       104.4       103.5          102.8
                ago.-19  104.1       104.2    103.7    103.7       104.5       105.4          104.5
                sep.-19  104.1       104.3    103.8    104.3       104.5       105.0          109.1

-Conforme a esto calcularemos los errores de cada uno de los modelos entre el valor observado de cada mes y la prediccion hecha de los modelos en el intervalo de 8 meses, reflejados en la tabla anterior.

-De todos las medidas que apareceran para determinar los errores de cada modelo, escogeremos RMSE por ser una de las mejores para determinar que pronostico, cuanto mas bajo salga el numero de dicha medida mejor pronostico sera.

-para poder comparar nuestros datos sera necesario almacenarlo en un objeto de serie de tiempo. De forma que:

real <- c(103.3, 103.7, 104.7, 104.9, 104.8, 104.2, 104.1, 104.1 )
data_real <- ts(real, frequency=12,start=c(2019,2))
data_real

##        Feb   Mar   Apr   May   Jun   Jul   Aug   Sep
## 2019 103.3 103.7 104.7 104.9 104.8 104.2 104.1 104.1

-Ahora haremos la comparacion de nuestros modelos con la data_real(valor de test en RMSE).

ARIMA:

accuracy(pron_real_ts,data_real)

##                 ME      RMSE       MAE      MPE      MAPE      ACF1 Theil's U
## Test set 0.1089892 0.2396797 0.2152161 0.104463 0.2065061 0.5099477 0.5125351

ETS:

accuracy(mod2,data_real)

##                      ME      RMSE       MAE       MPE      MAPE      MASE
## Training set 0.06916369 0.3084548 0.2472397 0.0681582 0.2449994 0.2247633
## Test set     0.45767849 0.4769796 0.4576785 0.4387901 0.4387901 0.4160714
##                   ACF1 Theil's U
## Training set 0.2998942        NA
## Test set     0.1159004  1.029046

HOLT-WINTERS:

accuracy(mod3,data_real)

##                 ME      RMSE       MAE      MPE      MAPE       ACF1 Theil's U
## Test set 0.3486004 0.4292439 0.4033664 0.334034 0.3866431 0.09569697 0.9323032

NNTEAR:

accuracy(mod4,data_real)

##                        ME      RMSE       MAE          MPE      MAPE      MASE
## Training set -0.000771296 0.5311222 0.3938052 -0.003679023 0.3893499 0.3580047
## Test set     -0.204681262 0.4554233 0.3941055 -0.197472445 0.3782473 0.3582778
##                   ACF1 Theil's U
## Training set 0.3646962        NA
## Test set     0.5877113  1.002075

-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:

m5 <- c(103.7, 103.4, 104.1, 102.6, 103.2, 103.5, 105.4, 105.0)
mod5c <- ts(m5, frequency=12,start=c(2019,2))
mod5c

##        Feb   Mar   Apr   May   Jun   Jul   Aug   Sep
## 2019 103.7 103.4 104.1 102.6 103.2 103.5 105.4 105.0

m6 <- c(104.0, 106.0, 106.6, 103.9, 103.2, 102.8, 104.5, 109.1)
mod6c <- ts(m6, frequency=12,start=c(2019,2))
mod6c

##        Feb   Mar   Apr   May   Jun   Jul   Aug   Sep
## 2019 104.0 106.0 106.6 103.9 103.2 102.8 104.5 109.1

-ELMAN(RNN):

accuracy(mod5c,data_real)

##              ME     RMSE    MAE       MPE      MAPE      ACF1 Theil's U
## Test set 0.3625 1.196349 1.0125 0.3441068 0.9692504 0.4660266  2.656639

JORDAN(RNN):

accuracy(mod6c,data_real)

##               ME     RMSE    MAE        MPE     MAPE      ACF1 Theil's U
## Test set -0.7875 2.238024 1.7875 -0.7592533 1.715147 0.2354782  4.999592

RESULTADOS:

           MODELO        RMSE (Test set)
           
           ARIMA           0.2396858
            ETS            0.4769796
           HOLT-W          0.4292113
           NNTEAR          0.4554233
         ELMAN(RNN)        1.196349
         JORDAN(RNN)       2.238024

-Por lo cual de todos los modelos sometidos a estudio el mejor modelo es ARIMA que tiene el menor RMSE 0.2396858