Roadmap for Time Series Analysis
PUJ/ Isabel Cristina García Arboleda: Carlos Alberto León Gil - Mauricio Pinzón Cortes
Mayo del 2023
Entregables y Avances
Avc1: Motivación
Objetivo: Identificar los desafíos actuales que implica la manipulación de observaciones generadas por medio de la variable tiempo.
Contexto
Actualmente los diversos productos financieros administrados en la compañía bancaria deben contar con diversos mecanismos de estimación que permitan mitigar los riesgos financieros asociados, por ello se plantea esta propuesta con el objetivo de identificar la relevancia de las diversas variables contempladas de negocio que permitan una estimación precisa y de alta confiabilidad. La entidad bancaria y la dirección de Business Intelligence son los principales afectados con la incertidumbre, la cual se asume debido a que no se cuentan con métodos de estimación a partir de los comportamientos históricos de los clientes, con la cual se permita proyectar los saldos en los siguientes periodos que permitan tomar decisiones estratégicas y de valor para los diversos productos administrados.
Descripción de los datos
Existen multitud de registros y atributos para el procesamiento dentro de las bases de datos transaccionales en la entidad bancaria, sin embargo, dentro de las consideraciones y delimitación de la propuesta se realizaron actividades primarias que se relacionan con el descubrimiento de las principales fuentes de información necesarias para el desarrollo de predicción de los saldos de productos de captación.
La mayoría de los tipos de valor de estos orígenes de datos se relacionan al detalle de los saldos y movimientos por cada cliente. Posterior a la actividad de investigación y descubrimiento de variables de interés, se realiza la selección de 89 variables numéricas con un total de 55.772 registros correspondientes a la transaccionalidad de 64 meses de clientes que poseen saldos mayor o igual a $100,000 COP al último día de cada mes, debido a las reglas de negocio que se establecen para la estimación de saldos a cierre. Si bien 64 de estas variables están relacionadas con el saldo final de cada cuenta al último día de cada mes, las 25 variables adicionales relacionan la descripción del cliente: edad, género, nivel educativo, estatus marital, personas a cargo, personas a cargo, tipo de vivienda, estrato, entre otros.
Para relacionar el comportamiento a nivel de negocio se definió por regla así como las autorizadas: código de subproducto, oficina, regional, zona, segmento comercial y declaración de renta.
Avc2: Estructura datos en TS
Objetivo Observar las tendencias y cambios estructurales en la serie a través de pruebas estadísticas para conocer la estructura subyacente de la serie.
Librerias
A continuación se enlistas las librerías usadas para el análisis de series temporales:
library(lubridate)
library(dplyr)
library(timeSeries)
library(forecast)
library(tseries)
library(ggplot2)
library(graphics)
require(magrittr)
require(dplyr)Visualización conjunto de datos
## # A tibble: 6 × 3
## Periodo_Original Periodo Saldo
## <chr> <dttm> <dbl>
## 1 Sum of SALDO_201801 2018-01-01 00:00:00 7683.
## 2 Sum of SALDO_201802 2018-02-01 00:00:00 7757.
## 3 Sum of SALDO_201803 2018-03-01 00:00:00 7845.
## 4 Sum of SALDO_201804 2018-04-01 00:00:00 7651.
## 5 Sum of SALDO_201805 2018-05-01 00:00:00 7699.
## 6 Sum of SALDO_201806 2018-06-01 00:00:00 8450.El set de datos origen contiene la información del saldo del periodo en estudio sin embargo, es importante seleccionar las dos variables con las cuales desarrollaremos el proceso de análisis de la serie de tiempo:
- Extracción de la fecha y saldo portafolio
head(Consolidado_Cuentas[,c(2,3)])## # A tibble: 6 × 2
## Periodo Saldo
## <dttm> <dbl>
## 1 2018-01-01 00:00:00 7683.
## 2 2018-02-01 00:00:00 7757.
## 3 2018-03-01 00:00:00 7845.
## 4 2018-04-01 00:00:00 7651.
## 5 2018-05-01 00:00:00 7699.
## 6 2018-06-01 00:00:00 8450.- Resumen del set de datos:
summary(Consolidado_Cuentas)## Periodo_Original Periodo Saldo
## Length:64 Min. :2018-01-01 00:00:00 Min. : 7651
## Class :character 1st Qu.:2019-04-23 12:00:00 1st Qu.: 9028
## Mode :character Median :2020-08-16 12:00:00 Median :11811
## Mean :2020-08-15 22:30:00 Mean :11002
## 3rd Qu.:2021-12-08 18:00:00 3rd Qu.:12866
## Max. :2023-04-01 00:00:00 Max. :13848Periodo: Tipo \(ymd\)
Se procede con el cambio de tipo de dato para la variable Periodo, ya que originalmente dentro del set de datos esta variable es de tipo \(character\) y mediante el paquete lubridate-package {lubridate} se realiza el cambio a \(ymd\) del campo Periodo.
## # A tibble: 6 × 2
## Periodo Saldo
## <date> <dbl>
## 1 2018-01-01 7683.
## 2 2018-02-01 7757.
## 3 2018-03-01 7845.
## 4 2018-04-01 7651.
## 5 2018-05-01 7699.
## 6 2018-06-01 8450.La función \(ts()\) convertirá un vector numérico el cual obedece al periodo de estudio de la serie temporal. Por otro lado, se identifica el periodo inicial y final que recibe la función como argumento en conjunto con la frecuencia que corresponde al número de observaciones antes de que se repita el patrón estacional en este caso la frecuencia es 12 (mensual).
- Identificación del periodo inicial y del periodo final de la serie de la serie de tiempo:
min(df$Periodo)## [1] "2018-01-01"max(df$Periodo)## [1] "2023-04-01"- Cambio de vector numérico a la serie de tiempo:
df.ts <- ts(data=df$Saldo, start=c(2018,1), end=c(2023,4), frequency = 12)
head(df.ts)## Jan Feb Mar Apr May Jun
## 2018 7683.489 7757.263 7845.047 7651.093 7699.410 8450.386modelo_red <- df.ts Visualización ts: Se identifica el comportamiento de la serie de tiempo de manera que se pueda visualizar el saldo portafolio para cada uno de los periodos
plot(df.ts, xlab="Años", ylab="Saldo")
title(main='Valores mensuales saldo portafolio')
Descriptiva: Un lag plot’ (gráfica de rezagos) se usa para evaluar cuando los valores de un conjunto de datos o una serie de tiempo es aleatoria. Si los datos son aleatorios, el lag plot no mostrará un comportamiento identificable, y si los datos no son aleatorios, el lag plot, mostrará un comportamiento identificable. Estos también, ayudan a identificar outliers.
lag.plot(df.ts, 12, do.lines = FALSE)
cycle(df.ts)## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2018 1 2 3 4 5 6 7 8 9 10 11 12
## 2019 1 2 3 4 5 6 7 8 9 10 11 12
## 2020 1 2 3 4 5 6 7 8 9 10 11 12
## 2021 1 2 3 4 5 6 7 8 9 10 11 12
## 2022 1 2 3 4 5 6 7 8 9 10 11 12
## 2023 1 2 3 4plot(aggregate(df.ts, FUN=mean))
boxplot(df.ts~cycle(df.ts))
Mediante el diagrama box plot logramos identificar que el mes 6 durante los últimos años presenta en promedio un crecimiento que se sostiene para finalmente cerrar diciembre con el promedio más alto en los saltos portafolio, esto puede ser el efecto del pago de las primas que se presenta en estos dos periodos de tiempo.
Promedio Móvil
Para realizar el promedio móvil se requiere el uso de la librería (forecast) el cual nos ayuda con el desarrollo de pronósticos basados en series de tiempo ya que este cuenta con la función de promedio móvil mediante \(ma\) Moving-average
Ma <- forecast::ma(df.ts, order = 12,centre = TRUE)
plot(df.ts)
lines(Ma, col="red")
Al calcular el promedio móvil se cumple el primer paso dentro de la lista de verificación para realizar la descomposición de las series de tiempo, en el cual se calcula la tendencia a largo plazo en los datos promediando entre las ventanas de datos para suavizar las fluctuaciones estacionales o irregularidades que existen dentro del comportamiento de los saldos de captación.
Avc3: Preprocesamiento y visualización
Objetivo: Manipular las bases de datos que involucran datos en series de tempo por medio de funciones para explicar el comportamiento regular en un intervalo de tiempo.
Descomposición Clasica
Si las fluctuaciones de los valores de la variable por encima y por debajo en la tendencia presentan un patrón irregular en el transcurso del tiempo, el modelo de agregación será aditivo. Si por el contrario la amplitud de dichas fluctuaciones crece o decrece con la tendencia el esquema o modelo es multiplicativo.
El procedimiento descomposición estacional indica dos métodos para modelar los factores estacionales
Estacionalidad Aditiva: Este método descomposición selecciona las correcciones estacionales y se agregan a la serie temporal corregida estacionalmente para obtener los valores observados. Estas correcciones tienen como objetivo eliminar de la serie el efecto estacional, todo esto con el propósito de estudiar aquellos comportamientos que quizás estén ocultos dentro de la serie de tiempo. Es de recordar que las observaciones son variación estacional obtendrán un componente estacional de 0.
seasonal_add <- df.ts - Ma
plot(seasonal_add)
Estacionalidad Multiplicativa: Este método selecciona el componente estacional con el cual se multiplica la serie corregida estacionalmente para dar lugar a la serie original. Las observaciones sin variación estacional tendrán un componente estacional de 1.
seasonal_multi <- df.ts/Ma
plot(seasonal_multi)
fit_add <- decompose(df.ts, type = c("additive", "multiplicative"))
plot(fit_add)
Como se puede observar mediante la gráfica de descomposición aditiva de la serie, en el panel superior encontramos los datos observados de nuestro comportamiento que hemos tratado todo el tiempo. En el siguiente panel encontramos la tendencia, el cual hemos calculado usando el promedio móvil (Moving-average smoothing) donde se aprecia un comportamiento similar a los datos observados.
En el tercer panel se aprecia la señal estacional y este comportamiento no tiene las fluctuaciones irregulares involucradas. Finalmente, en nuestro cuarto y último panel se identifica la señal aleatoria, es decir las fluctuaciones regulares en los datos que no se pueden explicar de nuestro set de datos estudiados.
fit_stl <- stl(df.ts, s.window = 12)
plot(fit_stl)
summary(fit_stl)## Call:
## stl(x = df.ts, s.window = 12)
##
## Time.series components:
## seasonal trend remainder
## Min. :-264.0189 Min. : 7562.703 Min. :-547.5591
## 1st Qu.:-112.5608 1st Qu.: 9048.984 1st Qu.: -65.6834
## Median : -37.4817 Median :11667.914 Median : 54.8074
## Mean : -2.2379 Mean :10974.168 Mean : 30.3882
## 3rd Qu.: 29.0804 3rd Qu.:12820.355 3rd Qu.: 181.4898
## Max. : 572.5543 Max. :13330.910 Max. : 397.0427
## IQR:
## STL.seasonal STL.trend STL.remainder data
## 141.6 3771.4 247.2 3837.8
## % 3.7 98.3 6.4 100.0
##
## Weights: all == 1
##
## Other components: List of 5
## $ win : Named num [1:3] 12 21 13
## $ deg : Named int [1:3] 0 1 1
## $ jump : Named num [1:3] 2 3 2
## $ inner: int 2
## $ outer: int 0Dickey–Fuller Test
Verificamos que sea ESTACIONARIA (Dicker-Fuller ‘adf.test’), si no lo es, entonces no podemos pronosticar más adelante usando modelos estacionarios.
test_ts <- df.ts
adf.test(test_ts)## Warning in adf.test(test_ts): p-value greater than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: test_ts
## Dickey-Fuller = 0.15694, Lag order = 3, p-value = 0.99
## alternative hypothesis: stationaryndiffs(df.ts)## [1] 1A partir del test Dickey Fuller obtenemos el \(p-value\) de 0.99, el cual nos indica que nuestro ts NO es estacionaria. Pensando en un nivel de significancia del 0.05 de nuevo, verifiquemos el número de veces que debemos diferenciar para tener la estacionariedad
diff_sea <- diff(df.ts, lag=1)
plot(diff_sea)
adf.test(diff_sea)## Warning in adf.test(diff_sea): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: diff_sea
## Dickey-Fuller = -4.1483, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationaryGraficos ACF y PACF
Como se aprecia con la segunda prueba de Dickey-Fuller test la serie temporal ya se encuentra estacionada mediante la diferenciación, con el cual el siguiente para ajustar un modelo ARIMA es determinar si se necesita términos AR o MA para corregir cualquier autocorrelación que permanezca en la serie diferenciada.
acf plot
ACF_1 <- acf(diff_sea)
plot(ACF_1)
pacf plot
PACF_1 <- pacf(diff_sea)
plot(PACF_1)
fit_stl_diff <- stl(diff_sea, s.window = 12)
plot(fit_stl_diff)
Avc4: Introducción a los modelos de pronóstico
Objetivo: Conocer los modelos básicos de series de tiempo usando las técnicas apropiadas y comunes para ejemplificar las relaciones entre este tipo de datos.
Modelo Holt - Winters
modelo_HW <- HoltWinters((diff_sea), seasonal = "additive")
plot(modelo_HW, main="Ajuste con Holt-Winters", xlab="Año", ylab="Saldo")
plot(fitted(modelo_HW), main='Descomposición', xlab='Año', ylab='Saldo')
pred <- predict(modelo_HW, 4, prediction.interval = TRUE)
pred## fit upr lwr
## May 2023 -227.6254 236.78426 -692.03503
## Jun 2023 557.1278 1027.16783 87.08771
## Jul 2023 -527.8380 -52.23415 -1003.44179
## Aug 2023 -516.7362 -35.63295 -997.83945plot(modelo_HW, pred)
fitted(modelo_HW)## xhat level trend season
## Feb 2019 113.703454 126.899158 -8.381415 -4.814289
## Mar 2019 46.029009 108.580406 -8.381415 -54.169982
## Apr 2019 -146.145628 90.246306 -8.381415 -228.010519
## May 2019 9.136303 73.248193 -8.381415 -55.730475
## Jun 2019 744.541358 59.832799 -8.381415 693.089974
## Jul 2019 -478.414894 49.582615 -8.381415 -519.616094
## Aug 2019 5.245029 40.832372 -8.381415 -27.205928
## Sep 2019 -2.223926 13.877852 -8.381415 -7.720364
## Oct 2019 7.968363 -17.236393 -8.381415 33.586171
## Nov 2019 -187.888376 -4.321718 -8.381415 -175.185243
## Dec 2019 393.179867 25.089144 -8.381415 376.472138
## Jan 2020 -5.371159 33.705645 -8.381415 -30.695389
## Feb 2020 -11.724963 26.070141 -8.381415 -29.413688
## Mar 2020 -60.184198 27.004593 -8.381415 -78.807376
## Apr 2020 -149.450779 108.271377 -8.381415 -249.340741
## May 2020 140.583718 217.156982 -8.381415 -68.191849
## Jun 2020 921.688783 241.606269 -8.381415 688.463928
## Jul 2020 -323.983848 204.926678 -8.381415 -520.529111
## Aug 2020 163.716667 245.280838 -8.381415 -73.182756
## Sep 2020 115.944445 188.320251 -8.381415 -63.994392
## Oct 2020 266.297642 188.375440 -8.381415 86.303617
## Nov 2020 36.389561 126.403254 -8.381415 -81.632278
## Dec 2020 575.245229 165.076983 -8.381415 418.549661
## Jan 2021 82.687398 119.917736 -8.381415 -28.848923
## Feb 2021 111.642258 126.376379 -8.381415 -6.352706
## Mar 2021 233.611750 98.880733 -8.381415 143.112432
## Apr 2021 98.333993 65.767280 -8.381415 40.948128
## May 2021 73.886289 69.188727 -8.381415 13.078977
## Jun 2021 664.942160 54.910420 -8.381415 618.413155
## Jul 2021 -351.214696 57.053259 -8.381415 -399.886540
## Aug 2021 -114.041147 87.778427 -8.381415 -193.438159
## Sep 2021 -25.934961 25.556440 -8.381415 -43.109987
## Oct 2021 -50.962472 3.776702 -8.381415 -46.357759
## Nov 2021 28.069771 1.600716 -8.381415 34.850470
## Dec 2021 327.465655 8.339166 -8.381415 327.507904
## Jan 2022 52.015921 52.510408 -8.381415 7.886928
## Feb 2022 -59.921431 2.129051 -8.381415 -53.669066
## Mar 2022 110.104094 36.596050 -8.381415 81.889459
## Apr 2022 32.636876 -29.147256 -8.381415 70.165547
## May 2022 -84.078183 -74.178272 -8.381415 -1.518496
## Jun 2022 533.850547 -102.233476 -8.381415 644.465438
## Jul 2022 -382.306221 -70.844731 -8.381415 -303.080074
## Aug 2022 -426.682682 -91.583362 -8.381415 -326.717905
## Sep 2022 -179.560846 -94.902540 -8.381415 -76.276891
## Oct 2022 -185.224843 -145.846908 -8.381415 -30.996520
## Nov 2022 -83.669083 -147.566637 -8.381415 72.278969
## Dec 2022 286.845880 -162.372186 -8.381415 457.599481
## Jan 2023 -291.204855 -186.741530 -8.381415 -96.081910
## Feb 2023 -166.887818 -210.906526 -8.381415 52.400123
## Mar 2023 -222.221062 -153.732502 -8.381415 -60.107145
## Apr 2023 -236.619269 -207.679074 -8.381415 -20.558779Suavizamiento Exponencial:
ts_ses <- ses(diff_sea, h = 4, level = c(80,95))
plot(ts_ses)
se ajusta un modelo de suavizamiento exponencial simple (SES) a la serie de tiempo diff_sea utilizando la función ses(). Se establece el parámetro de suavización alpha en default y se pronostican los próximos 4 meses utilizando el modelo ajustado.
accuracy(ts_ses)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -32.38316 325.7546 254.2376 386.1958 414.1873 1.162516 -0.08701336Avc5: Modelos estacionarios en series de tiempo
Objetivo: Esquematizar los modelos convencionales de series de tiempo a través de la metodología Box-Jenkins para encontrar patrones que nos permitan predecir futuras observaciones.
Modelos estacionarios en series de tiempo
Aplicaremos la metodología Box-Jenkins para identificar modelos autoregresivos integrados de media móvil (ARIMA) para analizar y predecir valores futuros de serie de tiempo.
Dentro de los pasos a seguir, tenemos
Visualizar la serie.
Transformarla en estacionaria.
Graficar ACF - PACF, escoger los parámetros.
Construir el modelo.
Hacer predicción.
Visualizar la serie.
plot((df.ts), xlab="Años", ylab="log(Saldo)")
title(main='Valores mensuales saldo portafolio')
- Transformarla en estacionaria.
plot(log(df.ts), xlab="Años", ylab="log(Saldo)")
title(main='Valores mensuales saldo portafolio')
- Transformarla en estacionaria.
Procedemos a transformar nuestra serie, para que sea estacionaria, haciendo uso de la función adf.
dif.saldo.ts <- (diff(log(df.ts)))
adf.test(dif.saldo.ts)##
## Augmented Dickey-Fuller Test
##
## data: dif.saldo.ts
## Dickey-Fuller = -4.2092, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationaryAl verificar mediante el test (Dicker-Fuller) si la serie de tiempo es estacionaria se evidencia que el valor p es menor que 0.05 con lo cual concluimos que la serie es estacionaria y procedemos a realizar el pronostico correspondiente
- Graficar ACF - PACF, escoger los parámetros.
ACF<-acf(dif.saldo.ts)
PACF<-pacf(dif.saldo.ts)
Le preguntamos a R, a través de la función auto.arima, cuál sería el modelo de ajuste. En el caso nuestro el modelo es ARIMA
modelo <- auto.arima(diff_sea)
modelo## Series: diff_sea
## ARIMA(0,1,1)(0,1,1)[12]
##
## Coefficients:
## ma1 sma1
## -0.7715 -0.5434
## s.e. 0.1036 0.2016
##
## sigma^2 = 57646: log likelihood = -346.54
## AIC=699.08 AICc=699.61 BIC=704.82el modelo ARIMA(0,0,0)(0,1,1)[12] implica que no hay componentes autoregresivos ni de media móvil en el modelo original de la serie de tiempo. Sin embargo, se ha incluido un término de media móvil en la serie de tiempo una vez diferenciada, y se ha tenido en cuenta la estacionalidad con un período de 12.
Ahora le pedimos los punto de cambio de la serie a R. La función se llama cpt.mean y detecta un punto de cambio en la media
require(changepoint)## Loading required package: changepoint## Loading required package: zoo##
## Attaching package: 'zoo'## The following object is masked from 'package:timeSeries':
##
## time<-## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Successfully loaded changepoint package version 2.2.4
## See NEWS for details of changes.mval<-cpt.mean(dif.saldo.ts, method = "AMOC")
cpts(mval)## numeric(0)modelo<-auto.arima(diff(log(df.ts)))
pronostico<- forecast(modelo,12,level=95)
plot(pronostico,main="Pronóstico con auto.arima",xlab="Año",ylab="log(Saldo)")
grid()
accuracy(pronostico)## ME RMSE MAE MPE MAPE MASE
## Training set 0.0008665001 0.01937703 0.01288142 -1.086085 152.1489 0.6293577
## ACF1
## Training set 0.1395081residuales<-modelo$residuals
qqnorm(residuales)
qqline(residuales)
shapiro.test(residuales)##
## Shapiro-Wilk normality test
##
## data: residuales
## W = 0.89273, p-value = 4.892e-05Box.test(residuales,type = "Ljung-Box")##
## Box-Ljung test
##
## data: residuales
## X-squared = 1.2855, df = 1, p-value = 0.2569Avc6: Regresión en series de tiempo
Objetivo: Combinar métodos autorregresivos y de aprendizaje estadístico mediante el enfoque de la clasificación para analizar datos de series de tiempo.
df_prophet <- data.frame(ds = as.Date(df$Periodo), origin = "2018-01-01", y = as.numeric(df$Saldo), by = 'm')
# Ajuste al modelo Prophet
model_prophet <- prophet(df_prophet)## Disabling weekly seasonality. Run prophet with weekly.seasonality=TRUE to override this.## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.future <- make_future_dataframe(model_prophet, periods = 12, freq = "m")
prophet_model <- predict(model_prophet, future)
# Visualizar los pronósticos
plot(model_prophet, prophet_model, xlab = 'Año', ylab = 'Saldo Captación')
Los valores pornosticados para el siguiente año es:
prophet_model$yhat## [1] 7607.651 7735.884 7745.645 7838.938 7871.327 8698.285 8469.409
## [8] 8353.455 8297.786 8417.478 8523.528 8946.557 8849.058 8927.940
## [15] 9045.030 8948.897 8843.134 9666.274 9266.067 9178.808 9006.722
## [22] 9144.200 9092.349 9698.742 9728.014 10064.586 10074.601 10534.689
## [29] 10845.038 11676.096 11790.272 11616.549 11789.493 11739.364 12022.135
## [36] 12517.929 12625.669 12743.867 12587.944 12805.323 12917.072 13607.062
## [43] 13414.579 13129.785 13049.135 13015.198 13138.712 13551.006 13506.671
## [50] 13570.164 13521.450 13361.313 13148.449 13721.978 13247.848 12878.465
## [57] 12569.367 12443.807 12296.431 12518.920 12214.197 12116.490 12074.079
## [64] 11801.407 11524.799 12171.399 11600.347 11336.549 10987.923 10952.334
## [71] 10721.648 10937.148 10562.499 10495.148 10127.299 10183.464splits <- initial_time_split(model_compare, prop = 0.90)
# Modelo 1: auto_arima
model_arima <- arima_reg() %>%
set_engine(engine = "auto_arima") %>%
fit(Saldo ~ Periodo, data = training(splits))## frequency = 12 observations per 1 year# Modelo 2: ets
model_ets <- exp_smoothing() %>%
set_engine(engine = "ets") %>%
fit(Saldo ~ Periodo, data = training(splits))## frequency = 12 observations per 1 year# Model 3: prophet
model_prophet <- prophet_reg(
seasonality_daily = FALSE,
seasonality_weekly = FALSE,
seasonality_yearly = TRUE) %>%
set_engine(engine = "prophet") %>%
fit(Saldo ~ Periodo, data = training(splits))
models_tbl <- modeltime_table(
model_arima,
model_ets,
model_prophet
)
# Calibración del modelo en el set de pruebas
calibration_tbl <- models_tbl %>%
modeltime_calibrate(new_data = testing(splits))
calibration_tbl <- models_tbl %>%
modeltime_calibrate(new_data = testing(splits))
calibration_tbl %>%
modeltime_forecast(
new_data = testing(splits),
actual_data = model_compare
) %>%
plot_modeltime_forecast(
.legend_max_width = 10,
.title = "Testeo de la Serie de Datos"
)refit_tbl <- calibration_tbl %>%
modeltime_refit(data = model_compare)## frequency = 12 observations per 1 year
## frequency = 12 observations per 1 yearrefit_tbl %>%
modeltime_forecast(h = "1 years", actual_data = model_compare) %>%
plot_modeltime_forecast(
.legend_max_width = 25
)calibration_tbl %>%
modeltime_accuracy() %>%
table_modeltime_accuracy(.interactive = FALSE)| Accuracy Table | ||||||||
| .model_id | .model_desc | .type | mae | mape | mase | smape | rmse | rsq |
|---|---|---|---|---|---|---|---|---|
| 1 | ARIMA(0,1,0)(1,0,0)[12] | Test | 657.47 | 5.44 | 2.57 | 5.27 | 723.81 | 0.00 |
| 2 | ETS(M,A,A) | Test | 287.75 | 2.39 | 1.13 | 2.35 | 347.38 | 0.38 |
| 3 | PROPHET | Test | 946.00 | 7.83 | 3.70 | 7.49 | 1006.17 | 0.12 |
Avc7: Redes neuronales recurrentes
Objetivo: Predecir valores futuros de las series de tiempo usando algoritmos propios de la inteligencia artificial para el tratamiento de las series de tiempo.
ELMAN NEURAL NETWORKS: La red neuronal de Elman es una red neuronal recurrente. En una red neuronal recurrente, las neuronas se conectan de nuevo con otras neuronas, el flujo de información es multidireccional, por lo que la activación de las neuronas puede fluir en un bucle. Este tipo de red neuronal tiene un sentido del tiempo y la memoria de los estados de las redes anteriores, lo que le permite aprender secuencias que varían con el tiempo, realizar tareas de clasificación y desarrollar predicciones de estados futuros. Como resultado, las redes neuronales recurrentes se utilizan para la clasificación, el modelado de secuencias estocásticas y las tareas de memoria asociativa.
Ahora, podemos estandarizar las observaciones usando el método de la escala
y <- as.ts(scale(modelo_red)) Tenga en cuenta que el método as.ts garantiza que conservemos un objeto ts.
Dado que no tenemos otros atributos explicativos para estos datos, utilizaremos un enfoque de series temporales puro. La pregunta principal en esta metodología de modelado es ¿cuántos retrasos de la variable dependiente usar? Bueno, ya que tenemos observaciones mensuales durante varios años, usemos 12 retrasos en la frecuencia mensual. Una forma de hacer esto es usar la operación Lag en el paquete quantmod. Esto requiere que convirtiéramos y en un objeto de clase de zool:
y<-as.zoo(y)
# Operador Quantmod Lag:
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)Esto nos da 12 atributos para alimentar como entradas en nuestra red neuronal.
SaldoCap_Neural <- cbind(x1,x2,x3,x4,x5,x6,x7,x8,x9, x10,x11,x12)
SaldoCap_Neural <- cbind(y,SaldoCap_Neural)
SaldoCap_Neural_Jordan <- cbind(y,SaldoCap_Neural)
head(round(SaldoCap_Neural,2),13)## Series 1 Lag.1 Lag.2 Lag.3 Lag.4 Lag.5 Lag.6 Lag.7 Lag.8 Lag.9 Lag.10 Lag.11
## 1 -1.63 NA NA NA NA NA NA NA NA NA NA NA
## 2 -1.59 -1.63 NA NA NA NA NA NA NA NA NA NA
## 3 -1.55 -1.59 -1.63 NA NA NA NA NA NA NA NA NA
## 4 -1.65 -1.55 -1.59 -1.63 NA NA NA NA NA NA NA NA
## 5 -1.62 -1.65 -1.55 -1.59 -1.63 NA NA NA NA NA NA NA
## 6 -1.25 -1.62 -1.65 -1.55 -1.59 -1.63 NA NA NA NA NA NA
## 7 -1.32 -1.25 -1.62 -1.65 -1.55 -1.59 -1.63 NA NA NA NA NA
## 8 -1.27 -1.32 -1.25 -1.62 -1.65 -1.55 -1.59 -1.63 NA NA NA NA
## 9 -1.23 -1.27 -1.32 -1.25 -1.62 -1.65 -1.55 -1.59 -1.63 NA NA NA
## 10 -1.16 -1.23 -1.27 -1.32 -1.25 -1.62 -1.65 -1.55 -1.59 -1.63 NA NA
## 11 -1.20 -1.16 -1.23 -1.27 -1.32 -1.25 -1.62 -1.65 -1.55 -1.59 -1.63 NA
## 12 -0.97 -1.20 -1.16 -1.23 -1.27 -1.32 -1.25 -1.62 -1.65 -1.55 -1.59 -1.63
## 13 -0.95 -0.97 -1.20 -1.16 -1.23 -1.27 -1.32 -1.25 -1.62 -1.65 -1.55 -1.59
## Lag.12
## 1 NA
## 2 NA
## 3 NA
## 4 NA
## 5 NA
## 6 NA
## 7 NA
## 8 NA
## 9 NA
## 10 NA
## 11 NA
## 12 NA
## 13 -1.63Se eliminan los NA de nuestro set de datos
SaldoCap_Neural <- SaldoCap_Neural[-(1:12),]
head(round(SaldoCap_Neural,2),13)## Series 1 Lag.1 Lag.2 Lag.3 Lag.4 Lag.5 Lag.6 Lag.7 Lag.8 Lag.9 Lag.10 Lag.11
## 13 -0.95 -0.97 -1.20 -1.16 -1.23 -1.27 -1.32 -1.25 -1.62 -1.65 -1.55 -1.59
## 14 -0.92 -0.95 -0.97 -1.20 -1.16 -1.23 -1.27 -1.32 -1.25 -1.62 -1.65 -1.55
## 15 -0.93 -0.92 -0.95 -0.97 -1.20 -1.16 -1.23 -1.27 -1.32 -1.25 -1.62 -1.65
## 16 -1.03 -0.93 -0.92 -0.95 -0.97 -1.20 -1.16 -1.23 -1.27 -1.32 -1.25 -1.62
## 17 -1.04 -1.03 -0.93 -0.92 -0.95 -0.97 -1.20 -1.16 -1.23 -1.27 -1.32 -1.25
## 18 -0.68 -1.04 -1.03 -0.93 -0.92 -0.95 -0.97 -1.20 -1.16 -1.23 -1.27 -1.32
## 19 -0.92 -0.68 -1.04 -1.03 -0.93 -0.92 -0.95 -0.97 -1.20 -1.16 -1.23 -1.27
## 20 -0.97 -0.92 -0.68 -1.04 -1.03 -0.93 -0.92 -0.95 -0.97 -1.20 -1.16 -1.23
## 21 -1.04 -0.97 -0.92 -0.68 -1.04 -1.03 -0.93 -0.92 -0.95 -0.97 -1.20 -1.16
## 22 -0.97 -1.04 -0.97 -0.92 -0.68 -1.04 -1.03 -0.93 -0.92 -0.95 -0.97 -1.20
## 23 -0.95 -0.97 -1.04 -0.97 -0.92 -0.68 -1.04 -1.03 -0.93 -0.92 -0.95 -0.97
## 24 -0.70 -0.95 -0.97 -1.04 -0.97 -0.92 -0.68 -1.04 -1.03 -0.93 -0.92 -0.95
## 25 -0.70 -0.70 -0.95 -0.97 -1.04 -0.97 -0.92 -0.68 -1.04 -1.03 -0.93 -0.92
## Lag.12
## 13 -1.63
## 14 -1.59
## 15 -1.55
## 16 -1.65
## 17 -1.62
## 18 -1.25
## 19 -1.32
## 20 -1.27
## 21 -1.23
## 22 -1.16
## 23 -1.20
## 24 -0.97
## 25 -0.95Estamos listos para comenzar a crear nuestros entrenamientos y pruebas. Primero, contamos el número de filas y utilizamos el método set.seed para la reproducibilidad
n <- nrow(SaldoCap_Neural)
set.seed(465)Como comprobación rápida, vemos que quedan 52 observaciones para su uso en el análisis. Esto es como se esperaba porque eliminamos 12 filas de observaciones debido al retraso en los datos. Usemos 39 filas de datos para construir el modelo, y las 13 filas de datos que utilizamos para la muestra de prueba. Seleccionamos los datos de entrenamiento al azar sin reemplazo de la siguiente manera
n_train <- 39
train <- sample(1:n,n_train , FALSE)Almacenamos los tributos de covariables que contienen los valores retrasados en las entradas del objeto R y la variable de respuesta en las salidas del objeto
inputs <- SaldoCap_Neural[,2:13]
outputs <- SaldoCap_Neural[,1]Encajamos una red neuronal con dos capas ocultas, cada una de las cuales contiene un nodo. Estimos la tasa de aprendizaje en 0,1 y el número máximo de iteraciones en 1000
fit <- elman(
inputs[train],
outputs[train],
size=c(1,1),
learnFuncParams=c(0.1),
maxit =1000
)Dado el tamaño relativamente pequeño del conjunto de datos, el modelo converge bastante rápido. Vamos a trazar la función de error:
plotIterativeError(fit)
El error cae muy rápido, nivelándose en alrededor de 400 y 600 iteraciones. También podemos utilizar el método de resumen para obtener más detalles sobre la red:
summary(fit)## SNNS network definition file V1.4-3D
## generated at Tue Jun 20 19:51:10 2023
##
## network name : RSNNS_untitled
## source files :
## no. of units : 17
## no. of connections : 20
## no. of unit types : 0
## no. of site types : 0
##
##
## learning function : JE_BP
## update function : JE_Order
##
##
## unit default section :
##
## act | bias | st | subnet | layer | act func | out func
## ---------|----------|----|--------|-------|--------------|-------------
## 1.00000 | 0.00000 | i | 0 | 1 | Act_Logistic | Out_Identity
## ---------|----------|----|--------|-------|--------------|-------------
##
##
## unit definition section :
##
## no. | typeName | unitName | act | bias | st | position | act func | out func | sites
## ----|----------|----------|----------|----------|----|----------|--------------|----------|-------
## 1 | | inp1 | 0.39372 | -0.52729 | i | 1, 1, 0 | Act_Identity | |
## 2 | | inp2 | 0.64615 | -0.15885 | i | 1, 2, 0 | Act_Identity | |
## 3 | | inp3 | 0.52197 | 0.76206 | i | 1, 3, 0 | Act_Identity | |
## 4 | | inp4 | 0.71463 | 0.83769 | i | 1, 4, 0 | Act_Identity | |
## 5 | | inp5 | 0.62402 | 0.25598 | i | 1, 5, 0 | Act_Identity | |
## 6 | | inp6 | 0.68532 | -0.20044 | i | 1, 6, 0 | Act_Identity | |
## 7 | | inp7 | 0.75534 | 0.83203 | i | 1, 7, 0 | Act_Identity | |
## 8 | | inp8 | 0.97738 | -0.65314 | i | 1, 8, 0 | Act_Identity | |
## 9 | | inp9 | 1.17103 | -0.30452 | i | 1, 9, 0 | Act_Identity | |
## 10 | | inp10 | 1.39766 | -0.93924 | i | 1,10, 0 | Act_Identity | |
## 11 | | inp11 | 1.01039 | 0.15395 | i | 1,11, 0 | Act_Identity | |
## 12 | | inp12 | 1.11355 | 0.48451 | i | 1,12, 0 | Act_Identity | |
## 13 | | hid11 | 0.74768 | 0.04192 | h | 7, 1, 0 |||
## 14 | | hid21 | 0.24514 | -6.56865 | h | 13, 1, 0 |||
## 15 | | out1 | 0.42519 | -900.67780 | o | 19, 1, 0 | Act_Identity | |
## 16 | | con11 | 1.07063 | 0.50000 | sh | 4,14, 0 | Act_Identity | |
## 17 | | con21 | 0.35799 | 0.50000 | sh | 10,14, 0 | Act_Identity | |
## ----|----------|----------|----------|----------|----|----------|--------------|----------|-------
##
##
## connection definition section :
##
## target | site | source:weight
## -------|------|---------------------------------------------------------------------------------------------------------------------
## 13 | | 16:-0.48832, 12: 1.20447, 11:-0.93896, 10: 0.58387, 9:-0.97484, 8:-0.02429, 7:-1.38142, 6: 1.16272, 5: 1.04670,
## 4:-0.92993, 3:-0.42277, 2: 1.41123, 1: 2.77745
## 14 | | 17: 0.23949, 13: 7.16070
## 15 | | 14: 1.73447
## 16 | | 16: 0.30000, 13: 1.00000
## 17 | | 17: 0.30000, 14: 1.00000
## -------|------|---------------------------------------------------------------------------------------------------------------------La primera columna proporciona un recuento del número de nodos o neuronas; vemos que toda la red tiene un total de 17. La tercera columna describe el tipo de neurona. Vemos que hay 12 neuronas de entrada, 2 neuronas ocultas, 2 neuronas en la capa de contexto y 1 neurona de salida. Las columnas cuarta y quinta dan el valor de la función de activación y el sesgo para cada neurona. Por ejemplo, la primera neurona tiene un valor de función de activación de \(21.0\) y un sesgo de \(-0.52729\)
Creando la predicción ideal
pred<-predict(fit, inputs[-train])
plot(pred)
El coeficiente de correlación al cuadrado:
cor(outputs[-train], pred)^2## [,1]
## [1,] 0.9123165Jordan Neural Networks: Las redes neuronales de Jordan son similares a la red neuronal de Elman. La única diferencia es que las neuronas de contexto se alimentan de la capa de salida en lugar de la capa oculta.
SaldoCap_Neural_Jordan <- SaldoCap_Neural_Jordan[-(1:12),]
plot(SaldoCap_Neural_Jordan)
Tenga en cuenta que la Serie 1 es la variable de respuesta y, y Lag 1, Lag 2,…,Lag 12 son los atributos de entrada x1, x2,…,x12
n=nrow(SaldoCap_Neural_Jordan)
n## [1] 52set.seed(465)Para la muestra de entrenamiento, 39 observaciones se muestran al azar sin reemplazo de la siguiente manera:
n_train <- 39
train <- sample(1:n,n_train , FALSE)
inputs <- SaldoCap_Neural_Jordan[ ,2:13]
outputs <- SaldoCap_Neural_Jordan[,1]
fit <- jordan(
inputs[train],
outputs[train],
size=2,
learnFuncParams=c(0.01),
maxit =1000)
plotIterativeError(fit)
pred <- predict(fit, inputs[-train])
cor(outputs[-train], pred)^2## [,1]
## [1,] 0.9322005