library(dplyr)
library(readxl)
library(lubridate)
library(nortest)
library(plotly)
library(stats)
library(tseries)
library(datasets)
library(base)
library(quantmod)
library(TSstudio)
library(forecast)Introducción
A lo largo del contenido se desarrolla el ejercicio planteado para la unidad de Modelos de Series de Tiempo del curso de Visualización Científica de la Maestría en Estadística Aplicada y la Maestría en Analítica de Datos de la Unuversidad del Norte.
Librerías
Las librerias utilizadas para el desarrollo del ejercicio propuesto son las siguientes:
API Yahoo Finence
Se considera la serie de tiempo asociada con las acciones de Tecnoglass desde que comenzó a comercializarse hasta la fecha del día de hoy.
acciones <- "TGLS"
TGLS <- getSymbols(acciones, src = "yahoo", auto.assign = FALSE)
TGLS_df <- data.frame(Date = index(TGLS),
Open = as.numeric(coredata(Op(TGLS))),
High = as.numeric(coredata(Hi(TGLS))),
Low = as.numeric(coredata(Lo(TGLS))),
Close = as.numeric(coredata(Cl(TGLS))),
Volume = as.numeric(coredata(Vo(TGLS))))Para empezar, se realiza la extracción de datos de Yahoo Finance, inicialmente como un data.frame para graficar el correspondiente candlestick:
plot_ly(data = TGLS_df, type = "candlestick",
x = ~Date,
open = ~Open,
close = ~Close,
high = ~High,
low = ~Low,
increasing = list(line = list(color = "#000D61")),
decreasing = list(line = list(color = "#DF6401"))) %>%
layout(title = "Gráfico de Velas de acciones TGLS - Tecnoglass",
xaxis = list(title = "Año"),
yaxis = list(title = "Precio"))En éste, se visualiza la evolución de los precios de apertura y cierre para cada uno de los días bursátiles considerados. Más específicamente, las velas azules indican que para ese día en particular, el precio de cierre fue más alto que el precio de apertura, mientras que las velas naranjas indican lo contrario: que el precio de cierre fue más bajo que el de apertura.
TGLS<- data.frame(Date = index(TGLS), Close = as.numeric(coredata(TGLS[, "TGLS.Close"])))
TGLS <- ts(TGLS$Close, start = c(year(min(TGLS$Date)), month(min(TGLS$Date))),end = c(year(max(TGLS$Date)), month(max(TGLS$Date))), frequency = 250)
ts_info(TGLS) The TGLS series is a ts object with 1 variable and 3000 observations
Frequency: 250
Start time: 2012 5
End time: 2024 4 El siguiente paso consiste en transformar el formato de los datos extraídos al formato de serie temporal, ya que este será el más adecuado para el desarrollo de lo propuesto. Se confirma de manera manual el formato del conjunto de datos: corresponde a una serie de tiempo (ts), contiene una sola variable y 3000 observaciones. La serie, en este caso, abarca desde mayo de 2012 hasta abril de 2024, con una frecuencia diaria.
Para este caso de estudio y para efectos de predecir, se incluyó en la serie de tiempo únicamente el precio de cierre de la acción TGLS.
Si bien es cierto que los datos tienen una frecuencia “diaría”, se establece frequency = 250 teniendo en cuenta que no todos los días son bursatiles, por lo que a los 365 días del año se debe restar los fines de semana y días festivos del país en el que la empresa cotiza. Tecnoglass cotiza en la Bolsa de Valores de Colombia (BVC) bajo el símbolo TGLSC. Sin embargo, para este informe se tomarón las acciones que se negocian en el mercado bursátil estadounidense, específicamente en el NASDAQ, bajo el símbolo TGLS.
Tendencia del precio de cierre
ts_plot(TGLS,
title = "TGLS Prices",
Ytitle = "Price",
Xtitle = "Year",
col = '#000D61')De acuerdo con lo mencionado de manera preliminar, nuestra variable de interés es el precio de cierre de la acción. Por esto, se grafica la tendencia de esta variable para obtener un bosquejo previo de la forma de los datos.
summary(TGLS) Min. 1st Qu. Median Mean 3rd Qu. Max.
2.290 8.848 10.180 14.762 14.910 59.000 El resumen nos muestra una amplia variabilidad, con un mínimo de 2.29 y un máximo de 59.00. La mediana se ubica en 10.18, mientras que la media es ligeramente más alta, alcanzando los 14.76. El primer cuartil, que indica el valor por debajo del cual se encuentra el 25% de los datos, es de 8.85, y el tercer cuartil, que indica el valor por debajo del cual se encuentra el 75% de los datos, es de 14.9.
Distribución del precio histórico de cierre
par(pty = "m", mar = c(5, 5, 2, 1) + 0.9)
hist(TGLS, main = "Histgram for TGLS",
xlab = "Precio de cierre",
breaks = "Sturges",
probability = TRUE,
col = 'pink')
lines(density(TGLS), col ='#000D61')
A través de la gráfica y en contraste con las medidas descritas en el paso anterior, podemos decir que la distribución de los precios de cierre de las acciones de Tecnoglass es asimétrica. La mayor parte de los precios de cierre se concentra en valores más bajos, con algunos valores extremadamente altos que contribuyen a una mayor media y un tercer cuartil más elevado en comparación con la mediana. La diferencia entre la media y la mediana también indica que la distribución está sesgada hacia la derecha.
Estacionariedad de la Serie Temporal
Al observar los datos de TGLS, se concluye una incostancia en su media y varianza a lo largo del tiempo, lo cual indica que la serie no posee estacionariedad. Esta falta de estacionariedad puede dificultar la modelización y predicción de la serie, ya que la media y la varianza cambian con el tiempo. Por lo tanto, es necesario explorar métodos de transformación que puedan abordar esta falta de estacionariedad para obtener resultados más precisos en el análisis de la serie de tiempo de TGLS.
Diferenciación de series temporales
Como primera opción, se decide aplicar el método de diferenciación, comenzando por la diferencuación de primer orden, como se plantea a continuación:
Diferencia de primer orden
ts_plot(diff(TGLS, lag = 1),
title = "TGLS Series - First Differencing",
Xtitle = "Year",
Ytitle = "Differencing of TICKER",
col = '#000D61')A través del método gráfico se obtiene una idea preliminar del comportamiento de la serie de tiempo. Asimismo, se continua el análisis, con otros diagramas y gráfico antes de proceder con la prueba de hipótesis.
ts_decompose(TGLS)La descomposición de la serie de tiempo nos muestra, en primer lugar, el comportamiento obsevado en la serie. Además, nos inidica que hasta el 2020 el comportamiento del precio de cierre de las acciones de Tecnoglass tenía una tendencia relativamente constante, sin embargo, a partir del año mencionado, la tendencia cambió al un alza sostenda, que apenas parece comenza a decrecer en los últimos peridos analizados.
acf(diff(TGLS, lag =1),
lag.max =250,
ylim = c(-0.5, 1),
col = '#000D61')
Con la primera diferenciación se evidencia estacionariedad en la serie, mediante el método gráfico. Sin embargo, para inferir si la serie es estacionaria o no, es conveniente realizar una prueba estadística.
Prueba de Hipótesis
\(\mathit{H_0}\) : La serie de tiempo para el precio de cierre de TGLS tras la primera diferencición no es estacionaria.
TGLS_D <- diff(TGLS, 1)
adf_result <- adf.test(TGLS_D)
adf_result
Augmented Dickey-Fuller Test
data: TGLS_D
Dickey-Fuller = -13.271, Lag order = 14, p-value = 0.01
alternative hypothesis: stationaryA partir de la prueba de Dickey-Fuller (\(\ est = -17.869, \ p-value = 0.01\)) se infiere estacionariedad en la serie tras la realizar la primera diferenciación.
par(mfrow=c(1,2))
acf(diff(TGLS), lag.max = 250, col = '#000D61')
pacf(diff(TGLS), lag.max = 250, col = '#000D61')
Elección del Modelo ARIMA - Criterio AIC
Para la elección de un modelo ARIMA óptimo a la hora de predecir, haremos uso de dos técnicas facilitadas por el docente durante el desarrollo de las clases. En primer lugar, se ejecutará un loop (código facilitado por el docente) que será contrastado con la función auto.arima del paquete forecast.
TGLS_split <- ts_split(TGLS, sample.out = 28)
train <- TGLS_split$train
test <- TGLS_split$testPara elegir un modelo eficaz y preciso, es necesario dividir la serie temporal en conjuntos de entrenamiento (train) y prueba (test). Podemos utilizar los datos de entrenamiento para ajustar un modelo y luego evaluar la precisión utilizando los datos de prueba. Esto nos permite verificar la capacidad del modelo para predecir correctamente sobre datos no incluidos en el conjunto de análisis.
par(mfrow=c(1,2))
acf(train, lag.max = 252, col = '#DF6401')
pacf(train, lag.max = 252, col = '#000D61')
Modelo ARIMA - Criterio AIC
Como método principal, se hace uso del loop mencionado previamente, creando una función en la que se utilizan los datos de entrenamiento y un límite para los posibles valores de los parámetros p, q y d.
best_ARIMA <- function(ts_in, p_n, d_n, q_n) {
best_aic <- Inf
best_pdq <- NULL
best_PDQ <- NULL
fit <- NULL
for(p in 1:p_n) {
for(d in 1:d_n) {
for (q in 1:q_n) {
for(P in 1:p_n) {
for(D in 1:d_n) {
for (Q in 1:q_n) {
tryCatch({
fit <- arima(scale(ts_in),
order=c(p, d, q),
seasonal = list(order = c(P, D, Q), period = 250),
xreg=1:length(ts_in),
method="CSS-ML")
tmp_aic <- AIC(fit)
if (tmp_aic < best_aic) {
best_aic <- tmp_aic
best_pdq = c(p, d, q)
best_PDQ = c(P, D, Q)
}
}, error=function(e){})
}
}
}
}
}
}
return(list("best_aic" = best_aic, "best_pdq" = best_pdq, "best_PDQ" = best_PDQ))
}best_model <- NULL
if(file.exists("best_arima.rda")) {
best_model = readRDS("best_arima.rda")
} else {
best_model = best_ARIMA(train, 4, 1, 4)
saveRDS(best_model, file = "best_arima.rda")
}Una vez definida la función con sus argumentos, se procede a determinar el modelo óptimo para nuestro conjunto de datos, sin olvidar que de manera previa se logró establecer que la serie de tiempo solo requiere una diferenciación para transformarse en estacionaria, por lo que el parámetro d se limita a 1.
best_model$best_aic
[1] -7167.355
$best_pdq
[1] 3 1 3
$best_PDQ
[1] 1 1 1El modelo obtenido a través del bucle, con parámetros order = (3, 1, 3) y seasonal = (1, 1, 1), tiene una parte no estacional de orden 3, lo que significa que toma en cuenta la relación entre los valores pasados hasta tres periodos anteriores. También tiene una parte estacional de orden 1, que considera cómo los datos varían en un patrón que se repite cada cierto período. Además, tanto para la parte no estacional como para la estacional, se tienen en cuenta los errores pasados hasta tres periodos atrás.
Modelo ARIMA - auto.arima
Como se ha dicho, el método auto.arima nos proporcionará información de contraste con el anterior.
if(file.exists("auto_arima.rda")) {
auto_model = readRDS("auto_arima.rda")
} else {
auto_model = auto.arima(train,
max.order = 3,
D = 1,
d = 1,
stepwise = FALSE,
approximation = FALSE)
saveRDS(auto_model, file = "auto_arima.rda")
}Cabe mencionar que, para este paso también se tiendrá en cuenta la limitación de los parámetros descrita en el punto anterior.
auto_modelSeries: train
ARIMA(3,1,0)(0,1,0)[250]
Coefficients:
ar1 ar2 ar3
0.0004 -0.0670 0.0558
s.e. 0.0192 0.0192 0.0193
sigma^2 = 0.6424: log likelihood = -3260.92
AIC=6529.84 AICc=6529.85 BIC=6553.47Y obtenemos un modelo en el que el 3 y el 1 representan los componentes del modelo que capturan cómo las observaciones recientes influyen en las futuras, mientras que el 0 indica que no se utilizan términos de medias móviles en este caso. El [252] señala que hay una tendencia estacional que se repite cada 252 días. La varianza del modelo es de 0.64, y la verosimilitud logarítmica sugiere un buen ajuste del modelo a los datos observados.
Entendido que a pesar de los buenos resultados obtenidos con la función auto.arima, se opta por el modelo proporcionado por el loop. Esto se debe a que se reconoce que la función auto.arima realiza una búsqueda corta, lo que podría llevar a omitir resultados que podrían dar lugar a un mejor modelo. El modelo elegido se presenta a continuación:
if(file.exists("model_arima.rda")) {
TGLS_model = readRDS("model_arima.rda")
} else {
TGLS_model <- arima(train, order = c(3,1,3),
seasonal = list(order = c(1,1,1)))
saveRDS(TGLS_model, file = "model_arima.rda")
}
TGLS_model
Call:
arima(x = train, order = c(3, 1, 3), seasonal = list(order = c(1, 1, 1)))
Coefficients:
ar1 ar2 ar3 ma1 ma2 ma3 sar1 sma1
-0.0438 0.0714 0.9348 0.0486 -0.1089 -0.9119 -0.1320 -0.7509
s.e. 0.0261 0.0251 0.0244 0.0298 0.0284 0.0281 0.0292 0.0251
sigma^2 estimated as 0.3935: log likelihood = -2725, aic = 5467.99Este modelo nos dice cómo las observaciones pasadas influyen en las futuras. Por ejemplo, si una observación pasada fue alta, es probable que la siguiente también lo sea, y viceversa. Además, este modelo también considera patrones estacionales en tus datos, como cambios que ocurren en ciertos períodos del año.
Los valores como ar1, ar2 y ar3 indican cómo las observaciones pasadas afectan a las futuras, mientras que ma1, ma2 y ma3 muestran cómo los errores pasados influyen en los futuros. Por otro lado, los términos como sar1 y sma1 explican cómo los errores de un período en la temporada influyen en los errores de los períodos siguientes.
El valor de sigma^2 nos dice cuánto varían nuestras predicciones con respecto a los datos reales. En este caso, un valor más bajo sugiere que nuestras predicciones están bastante cerca de los datos reales.
En resumen, este modelo ARIMA toma en cuenta tanto las influencias pasadas como las estacionales para predecir tus datos, y parece ajustarse bien a tus observaciones.
Predicción con Rolling
Para medir la precisión del modelo realizaremos predicciónes haciendo uso de una tecnica denominada Rolling. Ésta se utiliza para evaluar cómo se comporta un modelo a medida que se desplaza en el tiempo y cómo varía su precisión a medida que se agregan nuevos datos.
Además, implica entrenar el modelo repetidamente a medida que se avanza en el tiempo, utilizando solo la información disponible hasta ese momento para hacer la predicción siguiente. Esto significa que en cada paso de tiempo, se utiliza una ventana de datos que se desliza a lo largo de la serie temporal, y se entrena un modelo en esta ventana para hacer una predicción para el siguiente punto en el tiempo.
arima_rolling <- function(history, test, model) {
predictions <- c()
for (t in seq_along(test)) {
model_fit <- Arima(history, model=model)
forecast <- forecast(model_fit, h=1)
yhat <- forecast$mean
predictions <- c(predictions, yhat)
obs <- test[t]
history <- c(history, obs)
cat('predicted=', yhat, ', expected=', obs, '\n')
}
return(predictions)
}Utilizaremos la técnica seleccionada para predecir los siguientes horizontes: 7, 14, 21 y 28 días. Para los primeros tres horizontes se desplazará la ventana temporal de los datos de entrenamiento para que estos puedan complementar lo que se busca predecir.
Horizonte de 7 días
TGLS_split <- ts_split(TGLS, sample.out = 7)
train7 <- TGLS_split$train
test7 <- TGLS_split$testpredict7 <- arima_rolling(train7, test7, TGLS_model)predicted= 52.38442 , expected= 52.46
predicted= 51.87087 , expected= 53.94
predicted= 54.06421 , expected= 55.44
predicted= 55.34537 , expected= 58.41
predicted= 58.1515 , expected= 59
predicted= 59.64842 , expected= 59
predicted= 58.71942 , expected= 57.67 predict7[1] 52.38442 51.87087 54.06421 55.34537 58.15150 59.64842 58.71942train_df <- data.frame(Date = time(train7), Value = as.numeric(train7))
test_df <- data.frame(Date = time(test7), Value = as.numeric(test7))
predict_df <- data.frame(Date = time(test7), Value = predict7)
plot_ly() %>%
add_lines(data = train_df[-c(1:100), ], x = ~Date, y = ~Value, name = "Train", line = list(color = '#000D61')) %>%
add_lines(data = test_df, x = ~Date, y = ~Value, name = "Test", line = list(color = '#00C5DF')) %>%
add_lines(data = predict_df, x = ~Date, y = ~Value, name = "Forecast", line = list(color = '#DF6401')) %>%
layout(title = "Predicción TGLS - 7 días",
xaxis = list(title = "Fecha"),
yaxis = list(title = "Precio"),
showlegend = TRUE)Horizonte de 14 días
TGLS_split <- ts_split(TGLS, sample.out = 14)
train14 <- TGLS_split$train
test14 <- TGLS_split$testpredict14 <- arima_rolling(train14, test14, TGLS_model)predicted= 52.62159 , expected= 53.47
predicted= 53.86006 , expected= 53.6
predicted= 53.05593 , expected= 52.68
predicted= 52.85102 , expected= 52.03
predicted= 52.22682 , expected= 51.55
predicted= 51.53387 , expected= 50.75
predicted= 50.52599 , expected= 52.58
predicted= 52.38442 , expected= 52.46
predicted= 51.87087 , expected= 53.94
predicted= 54.06421 , expected= 55.44
predicted= 55.34537 , expected= 58.41
predicted= 58.1515 , expected= 59
predicted= 59.64842 , expected= 59
predicted= 58.71942 , expected= 57.67 predict14 [1] 52.62159 53.86006 53.05593 52.85102 52.22682 51.53387 50.52599 52.38442
[9] 51.87087 54.06421 55.34537 58.15150 59.64842 58.71942train_df <- data.frame(Date = time(train14), Value = as.numeric(train14))
test_df <- data.frame(Date = time(test14), Value = as.numeric(test14))
predict_df <- data.frame(Date = time(test14), Value = predict14)
plot_ly() %>%
add_lines(data = train_df[-c(1:100), ], x = ~Date, y = ~Value, name = "Train", line = list(color = '#000D61')) %>%
add_lines(data = test_df, x = ~Date, y = ~Value, name = "Test", line = list(color = '#00C5DF')) %>%
add_lines(data = predict_df, x = ~Date, y = ~Value, name = "Forecast", line = list(color = '#DF6401')) %>%
layout(title = "Predicción TGLS - 14 días",
xaxis = list(title = "Fecha"),
yaxis = list(title = "Precio"),
showlegend = TRUE)Horizonte de 21 días
TGLS_split <- ts_split(TGLS, sample.out = 21)
train21 <- TGLS_split$train
test21 <- TGLS_split$testpredict21 <- arima_rolling(train21, test21, TGLS_model)predicted= 45.89161 , expected= 45.46
predicted= 45.25383 , expected= 45.8
predicted= 46.2096 , expected= 45.42
predicted= 45.17621 , expected= 47.55
predicted= 47.54776 , expected= 50.43
predicted= 50.73506 , expected= 52.9
predicted= 52.54033 , expected= 52.26
predicted= 52.62159 , expected= 53.47
predicted= 53.86006 , expected= 53.6
predicted= 53.05593 , expected= 52.68
predicted= 52.85102 , expected= 52.03
predicted= 52.22682 , expected= 51.55
predicted= 51.53387 , expected= 50.75
predicted= 50.52599 , expected= 52.58
predicted= 52.38442 , expected= 52.46
predicted= 51.87087 , expected= 53.94
predicted= 54.06421 , expected= 55.44
predicted= 55.34537 , expected= 58.41
predicted= 58.1515 , expected= 59
predicted= 59.64842 , expected= 59
predicted= 58.71942 , expected= 57.67 predict21 [1] 45.89161 45.25383 46.20960 45.17621 47.54776 50.73506 52.54033 52.62159
[9] 53.86006 53.05593 52.85102 52.22682 51.53387 50.52599 52.38442 51.87087
[17] 54.06421 55.34537 58.15150 59.64842 58.71942train_df <- data.frame(Date = time(train21), Value = as.numeric(train21))
test_df <- data.frame(Date = time(test21), Value = as.numeric(test21))
predict_df <- data.frame(Date = time(test21), Value = predict21)
plot_ly() %>%
add_lines(data = train_df[-c(1:100), ], x = ~Date, y = ~Value, name = "Train", line = list(color = '#000D61')) %>%
add_lines(data = test_df, x = ~Date, y = ~Value, name = "Test", line = list(color = '#00C5DF')) %>%
add_lines(data = predict_df, x = ~Date, y = ~Value, name = "Forecast", line = list(color = '#DF6401')) %>%
layout(title = "Predicción TGLS - 14 días",
xaxis = list(title = "Fecha"),
yaxis = list(title = "Precio"),
showlegend = TRUE)Horizonte de 28 días
predict <- arima_rolling(train, test, TGLS_model)predicted= 42.9634 , expected= 44.25
predicted= 44.52705 , expected= 44.53
predicted= 44.60683 , expected= 45.01
predicted= 44.75045 , expected= 44.88
predicted= 44.64644 , expected= 44.65
predicted= 44.69806 , expected= 45
predicted= 45.06136 , expected= 45.58
predicted= 45.89161 , expected= 45.46
predicted= 45.25383 , expected= 45.8
predicted= 46.2096 , expected= 45.42
predicted= 45.17621 , expected= 47.55
predicted= 47.54776 , expected= 50.43
predicted= 50.73506 , expected= 52.9
predicted= 52.54033 , expected= 52.26
predicted= 52.62159 , expected= 53.47
predicted= 53.86006 , expected= 53.6
predicted= 53.05593 , expected= 52.68
predicted= 52.85102 , expected= 52.03
predicted= 52.22682 , expected= 51.55
predicted= 51.53387 , expected= 50.75
predicted= 50.52599 , expected= 52.58
predicted= 52.38442 , expected= 52.46
predicted= 51.87087 , expected= 53.94
predicted= 54.06421 , expected= 55.44
predicted= 55.34537 , expected= 58.41
predicted= 58.1515 , expected= 59
predicted= 59.64842 , expected= 59
predicted= 58.71942 , expected= 57.67 predict [1] 42.96340 44.52705 44.60683 44.75045 44.64644 44.69806 45.06136 45.89161
[9] 45.25383 46.20960 45.17621 47.54776 50.73506 52.54033 52.62159 53.86006
[17] 53.05593 52.85102 52.22682 51.53387 50.52599 52.38442 51.87087 54.06421
[25] 55.34537 58.15150 59.64842 58.71942train_df <- data.frame(Date = time(train), Value = as.numeric(train))
test_df <- data.frame(Date = time(test), Value = as.numeric(test))
predict_df <- data.frame(Date = time(test), Value = predict)
plot_ly() %>%
add_lines(data = train_df[-c(1:100), ], x = ~Date, y = ~Value, name = "Train", line = list(color = '#000D61')) %>%
add_lines(data = test_df, x = ~Date, y = ~Value, name = "Test", line = list(color = '#00C5DF')) %>%
add_lines(data = predict_df, x = ~Date, y = ~Value, name = "Forecast", line = list(color = '#DF6401')) %>%
layout(title = "Predicción TGLS - 28 días",
xaxis = list(title = "Fecha"),
yaxis = list(title = "Precio"),
showlegend = TRUE)accuracy(predict, test) ME RMSE MAE MPE MAPE ACF1 Theil's U
Test set 0.5297331 1.291709 0.9666665 1.027514 1.868746 0.2944323 1.028862Las predicciones y el resultado de las gráficas obtenidas a partir del uso del Rolling parecen indicar que obtivimos un modelo adecuado, preciso y ótimo en general; el análisis de las métricas de evaluación de la precisión del modelo en el conjunto de prueba, por su parte, muestran un contraste interesante.
Estos resultados indican que el modelo tiene un sesgo positivo leve en las predicciones, con un RMSE y MAE moderadamente bajos en comparación con los valores de los datos. Sin embargo, el MAPE es relativamente alto, lo que sugiere que el modelo puede estar cometiendo errores proporcionales significativos en algunas predicciones. Además, el valor de ACF1 sugiere que las predicciones pueden estar correlacionadas, lo que puede indicar que el modelo podría mejorarse mediante la incorporación de información adicional o mediante ajustes en los parámetros del modelo.
Predicción con forecast
La función forecast es útil para realizar pronósticos en nuestra serie de tiempo, adicionalmente, calcula intervalos de confianza alrededor de estos pronósticos.
Utilizaremos los mismos horizontes temporales planteados previamente:
Horizonte de 7 días
TGLS_7pred <- forecast(TGLS_model, h = 7)
TGLS_7pred Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
2023.904 42.96340 42.15917 43.76763 41.73343 44.19336
2023.908 43.23430 42.09423 44.37437 41.49071 44.97789
2023.912 43.35963 41.97954 44.73972 41.24896 45.47030
2023.916 43.06622 41.47214 44.66030 40.62828 45.50416
2023.920 42.84630 41.06349 44.62911 40.11972 45.57287
2023.924 42.93228 40.98975 44.87482 39.96143 45.90313
2023.928 42.95702 40.85949 45.05455 39.74912 46.16492accuracy(TGLS_7pred, test) ME RMSE MAE MPE MAPE MASE
Training set 0.01163525 0.6001854 0.2959441 -0.002172877 1.897691 0.06260861
Test set 1.79155121 1.8422250 1.7915512 3.987804939 3.987805 0.37901251
ACF1 Theil's U
Training set 0.01115154 NA
Test set 0.42707849 5.130238El modelo parece ajustarse bien al conjunto de entrenamiento, con errores medios y RMSE relativamente bajos. Sin embargo, en el conjunto de prueba, se observa un sesgo positivo en las predicciones, lo que indica que el modelo tiende a sobreestimar los valores. Además, el RMSE en el conjunto de prueba es bastante alto en comparación con el conjunto de entrenamiento, lo que sugiere que el modelo puede estar sobreajustado. El valor de Theil’s U también indica que las predicciones están significativamente alejadas de los valores reales. Esto sugiere que el modelo podría no generalizar bien a datos nuevos y desconocidos.
plot_forecast(TGLS_7pred,
title = "Close price TGLS - Forecast",
Ytitle = "Point Forecast",
Xtitle = "Year")Horizonte de 14 días
TGLS_14pred <- forecast(TGLS_model, h = 14)
TGLS_14pred Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
2023.904 42.96340 42.15917 43.76763 41.73343 44.19336
2023.908 43.23430 42.09423 44.37437 41.49071 44.97789
2023.912 43.35963 41.97954 44.73972 41.24896 45.47030
2023.916 43.06622 41.47214 44.66030 40.62828 45.50416
2023.920 42.84630 41.06349 44.62911 40.11972 45.57287
2023.924 42.93228 40.98975 44.87482 39.96143 45.90313
2023.928 42.95702 40.85949 45.05455 39.74912 46.16492
2023.932 43.29445 41.05342 45.53548 39.86709 46.72181
2023.936 43.13458 40.76645 45.50271 39.51284 46.75632
2023.940 43.47908 40.98427 45.97389 39.66360 47.29457
2023.944 43.29974 40.68578 45.91371 39.30203 47.29745
2023.948 43.27985 40.55773 46.00197 39.11672 47.44298
2023.952 43.62279 40.79153 46.45405 39.29276 47.95282
2023.956 43.36378 40.42920 46.29836 38.87573 47.85183accuracy(TGLS_14pred, test) ME RMSE MAE MPE MAPE MASE
Training set 0.01163525 0.6001854 0.2959441 -0.002172877 1.897691 0.06260861
Test set 3.49189844 4.4216710 3.4918984 7.174882097 7.174882 0.73873033
ACF1 Theil's U
Training set 0.01115154 NA
Test set 0.73790427 3.534842Estos resultados muestran que el modelo tiene un desempeño relativamente bueno en el conjunto de entrenamiento, con errores medios y RMSE bajos. Sin embargo, en el conjunto de prueba, se observa un sesgo positivo en las predicciones, indicando que tiende a sobreestimar los valores. El RMSE en el conjunto de prueba es bastante alto en comparación con el conjunto de entrenamiento, lo que sugiere que el modelo puede estar sobreajustado. Además, el valor de Theil's U indica que las predicciones están significativamente alejadas de los valores reales, lo que sugiere que el modelo puede no generalizar bien a nuevos datos. Sería importante considerar ajustes adicionales en el modelo o explorar otras técnicas de modelado para mejorar las predicciones.
plot_forecast(TGLS_14pred,
title = "Close price TGLS - Forecast",
Ytitle = "Point forecast",
Xtitle = "Year")Horizonte de 21 días
TGLS_21pred <- forecast(TGLS_model, h = 21)
TGLS_21pred Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
2023.904 42.96340 42.15917 43.76763 41.73343 44.19336
2023.908 43.23430 42.09423 44.37437 41.49071 44.97789
2023.912 43.35963 41.97954 44.73972 41.24896 45.47030
2023.916 43.06622 41.47214 44.66030 40.62828 45.50416
2023.920 42.84630 41.06349 44.62911 40.11972 45.57287
2023.924 42.93228 40.98975 44.87482 39.96143 45.90313
2023.928 42.95702 40.85949 45.05455 39.74912 46.16492
2023.932 43.29445 41.05342 45.53548 39.86709 46.72181
2023.936 43.13458 40.76645 45.50271 39.51284 46.75632
2023.940 43.47908 40.98427 45.97389 39.66360 47.29457
2023.944 43.29974 40.68578 45.91371 39.30203 47.29745
2023.948 43.27985 40.55773 46.00197 39.11672 47.44298
2023.952 43.62279 40.79153 46.45405 39.29276 47.95282
2023.956 43.36378 40.42920 46.29836 38.87573 47.85183
2023.960 43.72364 40.69367 46.75361 39.08971 48.35757
2023.964 44.09219 40.96538 47.21900 39.31015 48.87423
2023.968 43.68252 40.46375 46.90129 38.75984 48.60521
2023.972 43.81710 40.51230 47.12190 38.76285 48.87135
2023.976 44.00079 40.60841 47.39318 38.81258 49.18900
2023.980 44.08933 40.61364 47.56501 38.77372 49.40493
2023.984 43.82291 40.26844 47.37738 38.38682 49.25900accuracy(TGLS_21pred, test) ME RMSE MAE MPE MAPE MASE
Training set 0.01163525 0.6001854 0.2959441 -0.002172877 1.897691 0.06260861
Test set 5.15800458 6.1161554 5.1580046 10.176316778 10.176317 1.09120425
ACF1 Theil's U
Training set 0.01115154 NA
Test set 0.87844680 5.076679En este caso, lo resultados indican que el modelo tiene un buen desempeño en el conjunto de entrenamiento, con errores medios y RMSE bajos. Sin embargo, en el conjunto de prueba, se observa un sesgo positivo en las predicciones, lo que sugiere que el modelo tiende a sobreestimar los valores. El RMSE en el conjunto de prueba es bastante alto en comparación con el conjunto de entrenamiento, lo que indica que el modelo puede estar sobreajustado. Además, el valor de Theil's U sugiere que las predicciones están significativamente alejadas de los valores reales, lo que sugiere que el modelo puede no generalizar bien a nuevos datos. Sería importante considerar ajustes adicionales en el modelo o explorar otras técnicas de modelado
plot_forecast(TGLS_21pred,
title = "Close price TGLS - Forecast",
Ytitle = "Point forecast",
Xtitle = "Year")Horizonte de 28 días
TGLS_28pred <- forecast(TGLS_model, h = 28)
TGLS_28pred Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
2023.904 42.96340 42.15917 43.76763 41.73343 44.19336
2023.908 43.23430 42.09423 44.37437 41.49071 44.97789
2023.912 43.35963 41.97954 44.73972 41.24896 45.47030
2023.916 43.06622 41.47214 44.66030 40.62828 45.50416
2023.920 42.84630 41.06349 44.62911 40.11972 45.57287
2023.924 42.93228 40.98975 44.87482 39.96143 45.90313
2023.928 42.95702 40.85949 45.05455 39.74912 46.16492
2023.932 43.29445 41.05342 45.53548 39.86709 46.72181
2023.936 43.13458 40.76645 45.50271 39.51284 46.75632
2023.940 43.47908 40.98427 45.97389 39.66360 47.29457
2023.944 43.29974 40.68578 45.91371 39.30203 47.29745
2023.948 43.27985 40.55773 46.00197 39.11672 47.44298
2023.952 43.62279 40.79153 46.45405 39.29276 47.95282
2023.956 43.36378 40.42920 46.29836 38.87573 47.85183
2023.960 43.72364 40.69367 46.75361 39.08971 48.35757
2023.964 44.09219 40.96538 47.21900 39.31015 48.87423
2023.968 43.68252 40.46375 46.90129 38.75984 48.60521
2023.972 43.81710 40.51230 47.12190 38.76285 48.87135
2023.976 44.00079 40.60841 47.39318 38.81258 49.18900
2023.980 44.08933 40.61364 47.56501 38.77372 49.40493
2023.984 43.82291 40.26844 47.37738 38.38682 49.25900
2023.988 43.60069 39.96594 47.23545 38.04182 49.15957
2023.992 43.20217 39.49100 46.91335 37.52642 48.87793
2023.996 43.22238 39.43825 47.00651 37.43506 49.00970
2024.000 43.19031 39.33187 47.04876 37.28933 49.09130
2024.004 43.08937 39.16012 47.01862 37.08010 49.09864
2024.008 43.71748 39.72010 47.71486 37.60402 49.83095
2024.012 43.47221 39.40551 47.53892 37.25272 49.69170accuracy(TGLS_28pred, test) ME RMSE MAE MPE MAPE MASE
Training set 0.01163525 0.6001854 0.2959441 -0.002172877 1.897691 0.06260861
Test set 7.16940952 8.5534713 7.1694095 13.431440163 13.431440 1.51672803
ACF1 Theil's U
Training set 0.01115154 NA
Test set 0.89496258 6.46208Estos resultados indican que el modelo tiene un buen desempeño en el conjunto de entrenamiento, con errores medios y RMSE bajos. Sin embargo, en el conjunto de prueba, se observa un sesgo positivo en las predicciones, lo que sugiere que el modelo tiende a sobreestimar los valores. El RMSE en el conjunto de prueba es bastante alto en comparación con el conjunto de entrenamiento, lo que indica que el modelo puede estar sobreajustado. Además, el valor de Theil's U sugiere que las predicciones están significativamente alejadas de los valores reales, lo que sugiere que el modelo puede no generalizar bien a nuevos datos. Sería importante considerar ajustes adicionales en el modelo o explorar otras técnicas de modelado para mejorar las predicciones.
plot_forecast(TGLS_28pred, title = "Close price TGLS - Forecast", Ytitle = "Point forecast", Xtitle = "Year")test_forecast(TGLS,
forecast.obj = TGLS_28pred,
test = test)A medida que aumentamos el horizonte de predicción a 7, 14, 21 y 28 días, notamos un aumento significativo en la imprecisión de las predicciones. Esto se refleja en varios aspectos de la evaluación del modelo. Primero, tanto el error medio absoluto (MAE) como el error cuadrático medio de la raíz (RMSE) tienden a aumentar con el tiempo. Además, el error porcentual absoluto medio (MAPE) también aumenta, lo que indica una mayor discrepancia entre las predicciones y los valores reales a medida que avanzamos en el tiempo. Esto sugiere que el modelo tiene dificultades para capturar con precisión las tendencias y los patrones a largo plazo en los datos de la serie temporal. Además, el sesgo en las predicciones, medido por el error porcentual medio (MPE), tiende a aumentar, lo que indica una tendencia del modelo a sobreestimar o subestimar los valores a medida que nos alejamos del período de entrenamiento.
En resumen, estos resultados muestran que la precisión de las predicciones disminuye a medida que aumenta el horizonte de predicción, subrayando la importancia de ajustar y validar los modelos de series temporales para períodos más largos con el fin de obtener resultados más precisos y confiables.
Al comparar las predicciones obtenidas mediante el método de Rolling con las obtenidas mediante el método forecast, se nota una diferencia en su precisión. El método de Rolling parece ajustarse mejor a los datos de prueba, mostrando predicciones más cercanas a la realidad, incluso en horizontes de predicción más cortos como 7 días. Por otro lado, las predicciones generadas por el método forecast muestran una discrepancia más marcada con los datos reales, incluso en horizontes de predicción cortos. Esto sugiere que el método Rolling puede ser más consistente y confiable para predecir a corto plazo, mientras que el método forecast puede necesitar ajustes adicionales para mejorar su precisión en este contexto específico de predicción de series temporales.
Análisis de residuales
checkresiduals(TGLS_model)
Ljung-Box test
data: Residuals from ARIMA(3,1,3)(1,1,1)[250]
Q* = 1303.5, df = 492, p-value < 2.2e-16
Model df: 8. Total lags used: 500El valor p ($\ p-value = $) muy pequeño sugiere que hay evidencia significativa para rechazar la hipótesis nula de no autocorrelación en los residuos del modelo. Esto indica que el modelo ARIMA(3,1,3)(1,1,1)[250] no captura completamente la estructura de autocorrelación en los datos, lo que sugiere que podría haber espacio para mejorar el modelo o explorar otros enfoques para modelar la serie temporal.
Elección de Modelo ARIMA - criterio BIC*
best_ARIMA_BIC <- function(ts_in, p_n, d_n, q_n) {
best_bic <- Inf
best_pdq <- NULL
best_PDQ <- NULL
fit <- NULL
for(p in 1:p_n) {
for(d in 1:d_n) {
for (q in 1:q_n) {
for(P in 1:p_n) {
for(D in 1:d_n) {
for (Q in 1:q_n) {
tryCatch({
fit <- arima(scale(ts_in),
order=c(p, d, q),
seasonal = list(order = c(P, D, Q), period = 250),
xreg=1:length(ts_in),
method="CSS-ML")
tmp_bic <- BIC(fit)
if (tmp_bic < best_bic) {
best_bic <- tmp_bic
best_pdq <- c(p, d, q)
best_PDQ <- c(P, D, Q)
}
}, error=function(e){})
}
}
}
}
}
}
return(list("best_bic" = best_bic, "best_pdq" = best_pdq, "best_PDQ" = best_PDQ))
}best_model_bic <- NULL
if(file.exists("best_arima_bic.rda")) {
best_model_bic = readRDS("best_arima_bic.rda")
} else {
best_model_bic = best_ARIMA_BIC(train, 3, 1, 3)
saveRDS(best_model_bic, file = "best_arima_bic.rda")
}best_model_bic $best_bic
[1] -7108.267
$best_pdq
[1] 3 1 3
$best_PDQ
[1] 1 1 1Elección de Modelo ARIMA - criterio HQIC*
best_ARIMA_HQIC <- function(ts_in, p_n, d_n, q_n) {
best_hqic <- Inf
best_pdq <- NULL
best_PDQ <- NULL
fit <- NULL
for(p in 1:p_n) {
for(d in 1:d_n) {
for (q in 1:q_n) {
for(P in 1:p_n) {
for(D in 1:d_n) {
for (Q in 1:q_n) {
tryCatch({
fit <- arima(scale(ts_in),
order=c(p, d, q),
seasonal = list(order = c(P, D, Q), period = 250),
xreg=1:length(ts_in),
method="CSS-ML")
tmp_hqic <- AIC(fit, k = log(length(ts_in)))
if (tmp_hqic < best_hqic) {
best_hqic <- tmp_hqic
best_pdq <- c(p, d, q)
best_PDQ <- c(P, D, Q)
}
}, error=function(e){})
}
}
}
}
}
}
return(list("best_hqic" = best_hqic, "best_pdq" = best_pdq, "best_PDQ" = best_PDQ))
}best_model_hqic <- NULL
if(file.exists("best_arima_hqic.rda")) {
best_model_hqic = readRDS("best_arima_hqic.rda")
} else {
best_model_hqic = best_ARIMA_HQIC(train, 4, 1, 4)
saveRDS(best_model_hqic, file = "best_arima_hqic.rda")
}best_model_hqic $best_hqic
[1] -7107.385
$best_pdq
[1] 3 1 3
$best_PDQ
[1] 1 1 1Bajo el criterio BIC y HQIC, llegamos a la misma conclusión sobre el mejor modelo obtenida a través del criterio AIC. Esto indica que estos criterios evaluaron de manera similar la capacidad de los modelos para ajustarse a los datos y equilibrar su complejidad. Tanto el AIC y el BIC como el HQIC consideran tanto el ajuste del modelo como su complejidad, prefiriendo modelos más simples cuando dos modelos tienen un rendimiento similar en términos de ajuste. En este caso, la coincidencia en la selección del mejor modelo entre criterios sugiere que la elección del modelo final es confiable.
Conclusiones
Por supuesto, aquí tienes una versión equilibrada con un lenguaje más técnico pero aún accesible:
En sintesis, al analizar la serie de tiempo inicialmente encontramos que no era estacionaria, por lo que aplicamos la primera diferenciación para hacerla más manejable y poder modelar los datos adecuadamente. Probamos dos enfoques para seleccionar el mejor modelo: uno exhaustivo mediante un loop y otro utilizando la función auto.arima. Aunque ambos métodos nos dieron resultados, optamos por el modelo obtenido a través del bucle debido a la preocupación de que auto.arima pudiera pasar por alto aspectos importantes de los datos, lo que podría afectar la eficacia del modelo resultante.
Luego, hicimos predicciones utilizando las técnicas de Rolling y forecast. Descubrimos que el enfoque de rolling superó a forecast en términos de precisión y capacidad para capturar la variabilidad de los datos de prueba. Sin embargo, al evaluar las métricas y los residuos del modelo, descubrimos que el modelo seleccionado no era óptimo, lo que sugiere que aún hay espacio para mejorar la capacidad predictiva.
Es importante destacar que los criterios de información AIC, BIC y HQIC proporcionaron los mismos resultados en la elección del modelo óptimo, lo que respalda la consistencia de la evaluación del modelo independientemente del criterio utilizado.
Por último, este análisis resalta la importancia de contrastar varias técnicas de modelado y selección de modelos, además de la necesidad de evaluar la eficacia y el rendimiento de los modelos seleccionados. Aunque se han identificado puntos por mejorar, este ejercicio ofrece una base para el conocimiento y aprenizaje de modelos y análisis de series de tiempo.