library(dbplyr)
library(quantmod)
library(kableExtra)
library(TSstudio)
library(tidyverse)
library(yfR)
library(lubridate)
library(tseries)
library(plotly)
library(forecast)Para este informe, se utilizarán las siguientes librerías
Paso 1: Importación de datos y visualización
Considere la serie de tiempo asociada con las acciones de Tecnoglass desde que comenzó a comercializarse hasta la fecha del día de hoy.
Se importan los datos desde Yahoo Finance con el paquete yfR:
ticker <- "TGLS"
tecno0 <- getSymbols(ticker, src = "yahoo", auto.assign = FALSE)
tecno <- data.frame(ref_date = index(tecno0),
price_open = as.numeric(coredata(Op(tecno0))),
price_high = as.numeric(coredata(Hi(tecno0))),
price_low = as.numeric(coredata(Lo(tecno0))),
price_close = as.numeric(coredata(Cl(tecno0))),
volume = as.numeric(coredata(Vo(tecno0))))
tecno1_ts <- ts(tecno$price_close,
start = c(year(min(tecno$ref_date))),
end = c(year(max(tecno$ref_date))),
frequency = 250)
head(tecno, 10)Nota: Es importante anotar que se toma una frecuencia de 252 días debido a que la bolsa de valores no abre los fines de semana, además de tener en cuenta los años bisiestos entre 2012 y 2023 (3 años bisiestos).
A continuación, se visualiza la serie de tiempo del precio de cierre de Tecnoglass desde 2012 hasta 2024
tecno1 <- tecno%>%
select(ref_date, price_close)
ts_plot(tecno1_ts,
title = "Monthly Close Price 2012-2024",
Ytitle = 'Close price in USD',
Xtitle = 'Year',
color = '#006633')Se puede observar que la serie temporal presenta una tendencia a la alza tanto aditiva como multiplicativa a patir del año 2020, en donde se pasó de un cierre de 10 dólares a un cierre de aproximandamente 46 dólares. Esto muestra que la serie está en un estado no estacionario.
Resumen del precio de cierre:
summary(tecno1_ts) Min. 1st Qu. Median Mean 3rd Qu. Max.
2.29 8.85 10.18 14.78 14.91 59.00 A continuación, se muestra la frecuencia de los precios de cierre a lo largo de los años
hist(tecno1_ts, main = "Histogram for Close Price", xlab = "Close price", breaks = "Sturges", probability = TRUE, col = '#9999ff')
lines(density(tecno1$price_close), col = 'red')
A través del histograma, se puede evidenciar la estabilidad que había en el precio de cierre que había alrededor de los 10 dólares. Esto concuerda con el comportamiento constante de las acciones mostradas en la visualización de la seriey el candlestick mostrados previamente.
Además, puede observar que no hay datos faltantes para el precio de cierre
sum(is.na(tecno1_ts))[1] 0Se presenta un gráfico de vela del precio de cierre:
fig <- tecno %>% plot_ly(x = ~ref_date, type="candlestick",
open = ~price_open, close = ~price_close,
high = ~price_high, low = ~price_low)
fig <- fig %>% layout(title = "Close Price Candlestic",
xaxis = list(title = 'Day'),
yaxis = list(title = 'Close price in USD'))
figPaso 2: Proceso de predicción de resultados
Repita TODOS los pasos indicados en esta sección para encontrar modelos ARIMA para predecir el precio de las acciones de Tecnoglass con los siguientes horizontes: 7, 14 días, 21 días, 28 días. Utilizar siempre predicciones usando rolling con ventana de predicción continua de un día. Luego, use una ténica que no contenga rolling.
Para empezar, se aplica una prueba estadística para verificar si la serie de tiempo es estacionaria o no. Esta prueba es la de Dickey-Fuller. En esta prueba, se toman las siguientes hipótesis:
\[H_0: \textrm{La serie no es estacionaria}\] \[H_1: \textrm{La serie es estacionaria}\]
adf1 <- adf.test(tecno1_ts)
adf1
Augmented Dickey-Fuller Test
data: tecno1_ts
Dickey-Fuller = -0.10704, Lag order = 14, p-value = 0.99
alternative hypothesis: stationaryacf(tecno1_ts, lag.max = 120, ylim = c(-0.5, 1))
Teniendo en cuenta que el p-valor es de 0.99 > 0.05, no se rechaza la hipótesis nula, por lo que se puede concluir con un 95% de confianza que la serie temporal del precio de cierre no es estacionaria. Además, observando que la gráfica de ACF baja desde la parte positiva con una, aunque leve, tendencia lineal, se puede confirmar el comportamiento no estacionario obtenido con la prueba de Dickey Fuller.
Diferenciación de la serie temporal
Se realizan ahora figuras de autocorrelación para confirmar que la serie de tiempo diferenciada es estacionaria, así como también para verificar cuál es el orden de integración necesario para llevar la serie de tiempo no estacionaria a una estacionaría.
par(mfrow=c(2,2))
acf(diff(tecno1_ts), lag.max = 252, main = 'First differentation')
pacf(diff(tecno1_ts), lag.max = 252, main = 'First differentation')
acf(diff(diff(tecno1_ts)), lag.max = 252, main = 'Second differentation')
pacf(diff(diff(tecno1_ts)), lag.max = 252, main = 'Second differentation')
De las gráficas de ACF y PACF se puede observar cómo en la segunda diferenciación se entra en la zona negativa con gran rapidez. Esto muestra cómo la serie se encuentra sobrediferenciada y, por lo tanto, es mejor optar por realizar una sola diferenciación.
Es importante resaltar que estas con conclusiones obtenidas por un análisis visual, y por eso, es importantre usar los criterios de AIC y BIC como mecanismo para decidir con certeza cual es el número óptimo de diferenciaciones que deben aplicar a la serie temporal.
Para comprobar la estacionareidad de la serie, se aplica el test de Dickey-Fuller.
adf.test(diff(tecno1_ts, 1))
Augmented Dickey-Fuller Test
data: diff(tecno1_ts, 1)
Dickey-Fuller = -13.209, Lag order = 14, p-value = 0.01
alternative hypothesis: stationaryTeniendo en cuenta que el p-valor es de 0.99 < 0.05, no se rechaza la hipótesis nula, por lo que se puede concluir con un 95% de confianza que la serie temporal del precio de cierre es estacionaria. Por lo tanto, es sufuciente aplicar diferenciación de primer orden para obtener una estabilidad en la serie de precio de cierre.
A continuación se pueden visualizar las series de tiempo con las diferenciaciones de cada orden:
ts_plot(diff(tecno1_ts, lag = 1),
title = "Close Price Series - First Differencing",
Xtitle = "Year",
Ytitle = "Differencing of Close Price in USD",
color = "#006633")Proceso de entrenamiento
Antes de aplicar los criterios de AIC o BIC es necesario hacer las particiones de training y testing de los datos usando la función ts_split.
ts_info(tecno1_ts) The tecno1_ts series is a ts object with 1 variable and 3001 observations
Frequency: 250
Start time: 2012 1
End time: 2024 1 A continuación se muestra la descomposición de las serie de tiempo aditivas:
ts_decompose(tecno1_ts)Se define la partición para los diferentes un horizonte de 28 días, siendo este el mayor tiempo de los 4 requeridos. Este mismo entrenamiento se utilizará para las diferentes predicciones:
Horizonte de 28 días
tecno_split_28 <- ts_split(tecno1_ts, sample.out = 28)
train28 <- tecno_split_28$train
test28 <- tecno_split_28$test
par(mfrow=c(1,2))
acf(train28, lag.max = 60)
pacf(train28, lag.max = 60)
De la gráfica de autocorrelación se puede observar cómo el conjunto de traning presenta un comportamiento lineal, por lo que no se tiene la condición de estacionareidad. Este problema se arreglará internamente en la generación del modelo. Sin embargo, este resultado es importante para determinar el parámetro de diferenciación de dicho modelo.
Se realiza el test de Dickey-Fuller para comprobar que a través de la diferenciación es posible conseguir la condición de estacionareidad.
adf.test(diff(train28, 1))
Augmented Dickey-Fuller Test
data: diff(train28, 1)
Dickey-Fuller = -13.506, Lag order = 14, p-value = 0.01
alternative hypothesis: stationaryAl obtener un p-valor menor a 0.05, se puede concluir que se rechaza la hipótesis de no estacionareidad, lo cual permite fijar un parámetro de diferenciación \(d = 1\) al momento de calcular el modelo ARIMA.
Criterio AIC
Ahora se que tiene un panorama claro del comportamiento del entramiento de la serie, se procede a calcular el modelo ARIMA óptimo para realizar las predicciones. Se empieza con el método adoptado por la función auto.arima, el cual minimiza el puntaje AIC, maximizando la máxima verosimilitud.
if(file.exists("auto_arima.rda")) {
auto_model = readRDS("auto_arima.rda")
} else {
auto_model = auto.arima(train28,
max.order = 3,
D = 1,
d = 1,
stepwise = FALSE,
approximation = FALSE)
saveRDS(auto_model, file = "auto_arima.rda")
}
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.47A través de esta función, se obtiene que el modelo ARIMA con el menor puntaje AIC es ARIMA(0,1,2)[0,1,0] con un puntaje de 6530.
Antes de implementar este modelo, se opta por probar una manera distinta de hallar el modelo la cual encuentra el mejor modelo ARIMA con base en la minimización del coeficiente de Akaike (AIC) para distintos rangos de los parámetros \(p,d,q\) del modelo, para obtener el mejor modelo en términos de bondad de ajuste. Esto evita los posibles problemas que genera la función auto.arima que puede pasar por alto mejores modelos.
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))
}
if(file.exists("best_arima.rda")) {
best_model = readRDS("best_arima.rda")
} else {
best_model = best_ARIMA(train28, 3, 1, 3)
saveRDS(best_model, file = "best_arima.rda")
}
best_model$best_aic
[1] -7167.355
$best_pdq
[1] 3 1 3
$best_PDQ
[1] 1 1 1A través de este ciclo, se obtuvo que el mejor modelo ARIMA es ARIMA(3,1,3)[1,1,1] con un puntaje AIC de -7167. En estos estudios, se ignora el signo del puntaje y simplemente se escoge el modelo con un menor puntaje. En este caso, el menor puntaje fue obtenido por el modelo generado por la función auto.arima.
Siendo así, se procede a recontruir el modelo con los mejores parámetros.
if(file.exists("TGLS_model.rda")) {
fit_model = readRDS("TGLS_model.rda")
} else {
fit_model = arima(train28, order = c(3,1,3),
seasonal = list(order = c(1,1,1)))
saveRDS(fit_model, file = "TGLS_model.rda")
}
fit_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.99Predicciones (Con Rolling)
Teniendo el modelo entrenado y con los parámetros óptimos, es momento de realizar las predicciones de la serie. En este apartado, se utiliza en particular el método de Rolling.
Primero, se crea la función con la que se ejecutará cada una de las predicciones:
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)
}Ahora, se hace la partición de training y testing para cada uno de los horizontes requeridos por el ejercicio:
Horizonte de 7 días:
tecno1_split <- ts_split(tecno1_ts, sample.out = 7)
train7 <- tecno1_split$train
test7 <- tecno1_split$testpredict7 <- arima_rolling(train7, test7, fit_model)predicted= 51.87087 , expected= 53.94
predicted= 54.06421 , expected= 55.44
predicted= 55.34537 , expected= 58.41
predicted= 58.15151 , expected= 59
predicted= 59.64842 , expected= 59
predicted= 58.71942 , expected= 57.67
predicted= 57.6202 , expected= 58.72 predict7[1] 51.87087 54.06421 55.34537 58.15151 59.64842 58.71942 57.62020train_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 = '#006633')) %>%
add_lines(data = test_df, x = ~Date, y = ~Value, name = "Test", line = list(color = 'red')) %>%
add_lines(data = predict_df, x = ~Date, y = ~Value, name = "Forecast", line = list(color = 'navy')) %>%
layout(title = "Tecnoglass 7 days predictions",
xaxis = list(title = "Date"),
yaxis = list(title = "Close Price in USD"),
showlegend = TRUE)accuracy(predict7, test7) ME RMSE MAE MPE MAPE ACF1 Theil's U
Test set 0.9657138 1.648257 1.450812 1.708097 2.542017 0.3394059 1.014421Horizonte de 14 días:
tecno1_split <- ts_split(tecno1_ts, sample.out = 14)
train14 <- tecno1_split$train
test14 <- tecno1_split$testpredict14 <- arima_rolling(train14, test14, fit_model)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.15151 , expected= 59
predicted= 59.64842 , expected= 59
predicted= 58.71942 , expected= 57.67
predicted= 57.6202 , expected= 58.72 predict14 [1] 53.86006 53.05593 52.85102 52.22682 51.53387 50.52599 52.38442 51.87087
[9] 54.06421 55.34537 58.15151 59.64842 58.71942 57.62020train_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 = '#006633')) %>%
add_lines(data = test_df, x = ~Date, y = ~Value, name = "Test", line = list(color = 'red')) %>%
add_lines(data = predict_df, x = ~Date, y = ~Value, name = "Forecast", line = list(color = 'navy')) %>%
layout(title = "Tecnoglass 14 days predictions",
xaxis = list(title = "Date"),
yaxis = list(title = "Close Price in USD"),
showlegend = TRUE)accuracy(predict14, test14) ME RMSE MAE MPE MAPE ACF1 Theil's U
Test set 0.4265642 1.341574 1.085926 0.740924 1.962779 0.2853145 1.074159Horizonte de 21 días:
tecno1_split <- ts_split(tecno1_ts, sample.out = 21)
train21 <- tecno1_split$train
test21 <- tecno1_split$testpredict21 <- arima_rolling(train21, test21, fit_model)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.15151 , expected= 59
predicted= 59.64842 , expected= 59
predicted= 58.71942 , expected= 57.67
predicted= 57.6202 , expected= 58.72 predict21 [1] 45.25383 46.20960 45.17621 47.54776 50.73506 52.54033 52.62159 53.86006
[9] 53.05593 52.85102 52.22682 51.53387 50.52599 52.38442 51.87087 54.06421
[17] 55.34537 58.15151 59.64842 58.71942 57.62020train_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 = '#006633')) %>%
add_lines(data = test_df, x = ~Date, y = ~Value, name = "Test", line = list(color = 'red')) %>%
add_lines(data = predict_df, x = ~Date, y = ~Value, name = "Forecast", line = list(color = 'navy')) %>%
layout(title = "Tecnoglass 21 days predictions",
xaxis = list(title = "Date"),
yaxis = list(title = "Close Price in USD"),
showlegend = TRUE)accuracy(predict21, test21) ME RMSE MAE MPE MAPE ACF1 Theil's U
Test set 0.653215 1.472716 1.194687 1.22273 2.253952 0.2796751 1.029434Horizonte de 28 días:
predict28 <- arima_rolling(train28, test28, fit_model)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.15151 , expected= 59
predicted= 59.64842 , expected= 59
predicted= 58.71942 , expected= 57.67
predicted= 57.6202 , expected= 58.72 predict28 [1] 44.52705 44.60683 44.75045 44.64644 44.69806 45.06136 45.89161 45.25383
[9] 46.20960 45.17621 47.54776 50.73506 52.54033 52.62159 53.86006 53.05593
[17] 52.85102 52.22682 51.53387 50.52599 52.38442 51.87087 54.06421 55.34537
[25] 58.15151 59.64842 58.71942 57.62020train_df <- data.frame(Date = time(train28), Value = as.numeric(train28))
test_df <- data.frame(Date = time(test28), Value = as.numeric(test28))
predict_df <- data.frame(Date = time(test28), Value = predict28)
plot_ly() %>%
add_lines(data = train_df[-c(1:100), ], x = ~Date, y = ~Value, name = "Train", line = list(color = '#006633')) %>%
add_lines(data = test_df, x = ~Date, y = ~Value, name = "Test", line = list(color = 'red')) %>%
add_lines(data = predict_df, x = ~Date, y = ~Value, name = "Forecast", line = list(color = 'navy')) %>%
layout(title = "Tecnoglass 28 days predictions",
xaxis = list(title = "Date"),
yaxis = list(title = "Close Price in USD"),
showlegend = TRUE)accuracy(predict28, test28) ME RMSE MAE MPE MAPE ACF1 Theil's U
Test set 0.5230617 1.285531 0.9599951 0.990563 1.831795 0.2833021 1.030445Al observar los resultados para los diferentes horizontes de predicción, notamos una mejora gradual en varias métricas de evaluación a medida que aumenta el horizonte de predicción. Para el horizonte de 7 días, la media del error (ME) es positiva pero pequeña (0.97), lo que indica un ligero sesgo hacia la sobreestimación. El error cuadrático medio (RMSE) y el error absoluto medio (MAE) son más altos, lo que sugiere una dispersión relativamente grande de los errores de predicción. Además, los errores porcentuales (MPE y MAPE) son moderados, indicando una precisión razonable en las predicciones.
Para el horizonte de 14 días, observamos una mejora en varias métricas de evaluación en comparación con el horizonte de 7 días. La ME sigue siendo positiva pero disminuye, lo que indica un sesgo menor en las predicciones. El RMSE y el MAE también disminuyen, lo que sugiere una reducción en la dispersión de los errores. Los errores porcentuales (MPE y MAPE) también muestran una disminución, indicando una mayor precisión en las predicciones en comparación con el horizonte de 7 días.
Para el horizonte de 21 días, continuamos viendo mejoras en las métricas de evaluación en comparación con los horizontes anteriores. La ME sigue siendo positiva, aunque ligeramente mayor que para el horizonte de 14 días. El RMSE y el MAE muestran una tendencia a la baja, lo que indica una disminución adicional en la dispersión de los errores. Los errores porcentuales (MPE y MAPE) también muestran una ligera disminución, lo que sugiere una mayor precisión en las predicciones en comparación con los horizontes anteriores.
Finalmente, para el horizonte de 28 días, observamos una continuación de las tendencias positivas en las métricas de evaluación. La ME sigue siendo positiva pero disminuye aún más en comparación con los horizontes anteriores. El RMSE y el MAE muestran una disminución adicional, lo que indica una reducción continua en la dispersión de los errores. Los errores porcentuales (MPE y MAPE) también disminuyen, lo que sugiere una mayor precisión en las predicciones en comparación con los horizontes anteriores.
En resumen, estos resultados sugieren que el modelo de predicción mejora gradualmente su desempeño a medida que se proyecta a horizontes más largos, lo que respalda la observación inicial de que las predicciones se acercan más al valor observado en el conjunto de prueba a medida que aumenta el horizonte de predicción. Sin embargo, para confirmar la confiabilidad de estos resultados, es crucial realizar un análisis de los residuos del modelo.
Predicciones (Sin Rolling)
Para esta parte, se realizan las predicciones sin utilizar el método de Rolling. En este caso, se toma la función forecast.
Horizonte de 7 días
tecno7_pred <- forecast(auto_model, h = 7)
tecno7_pred Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
2023.904 42.81117 41.78404 43.83830 41.24031 44.38203
2023.908 43.20833 41.75545 44.66122 40.98634 45.43033
2023.912 42.86101 41.12029 44.60174 40.19880 45.52322
2023.916 40.34931 38.33378 42.36484 37.26682 43.43179
2023.920 40.55486 38.29564 42.81408 37.09968 44.01003
2023.924 40.35932 37.88339 42.83524 36.57272 44.14592
2023.928 38.61063 35.93441 41.28686 34.51771 42.70356plot_forecast(tecno7_pred,
title = "7 days Close price - Forecast",
Ytitle = "Point Forecast",
Xtitle = "Year")accuracy(tecno7_pred, test28) ME RMSE MAE MPE MAPE MASE
Training set 0.0002181155 0.7664593 0.3852712 -0.01923693 2.450196 0.08150625
Test set 4.3864805159 4.9733242 4.3864805 9.52877864 9.528779 0.92798408
ACF1 Theil's U
Training set 0.001789564 NA
Test set 0.364122289 5.664439Horizonte de 14 días
tecno14_pred <- forecast(auto_model, h = 14)
tecno14_pred Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
2023.904 42.81117 41.78404 43.83830 41.24031 44.38203
2023.908 43.20833 41.75545 44.66122 40.98634 45.43033
2023.912 42.86101 41.12029 44.60174 40.19880 45.52322
2023.916 40.34931 38.33378 42.36484 37.26682 43.43179
2023.920 40.55486 38.29564 42.81408 37.09968 44.01003
2023.924 40.35932 37.88339 42.83524 36.57272 44.14592
2023.928 38.61063 35.93441 41.28686 34.51771 42.70356
2023.932 40.33076 37.46795 43.19357 35.95247 44.70905
2023.936 37.84036 34.80260 40.87812 33.19451 42.48621
2023.940 39.32043 36.11723 42.52362 34.42156 44.21929
2023.944 39.28046 35.91994 42.64099 34.14098 44.41994
2023.948 38.45043 34.93963 41.96123 33.08113 43.81974
2023.952 38.70043 35.04553 42.35533 33.11075 44.29012
2023.956 38.92044 35.12691 42.71397 33.11873 44.72215plot_forecast(tecno14_pred,
title = "14 days Close price - Forecast",
Ytitle = "Point Forecast",
Xtitle = "Year")accuracy(tecno14_pred, test28) ME RMSE MAE MPE MAPE MASE
Training set 0.0002181155 0.7664593 0.3852712 -0.01923693 2.450196 0.08150625
Test set 8.9451463760 10.2382297 8.9451464 17.60644629 17.606446 1.89239493
ACF1 Theil's U
Training set 0.001789564 NA
Test set 0.820640035 7.575408Horizonte de 21 días
tecno21_pred <- forecast(auto_model, h = 21)
tecno21_pred Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
2023.904 42.81117 41.78404 43.83830 41.24031 44.38203
2023.908 43.20833 41.75545 44.66122 40.98634 45.43033
2023.912 42.86101 41.12029 44.60174 40.19880 45.52322
2023.916 40.34931 38.33378 42.36484 37.26682 43.43179
2023.920 40.55486 38.29564 42.81408 37.09968 44.01003
2023.924 40.35932 37.88339 42.83524 36.57272 44.14592
2023.928 38.61063 35.93441 41.28686 34.51771 42.70356
2023.932 40.33076 37.46795 43.19357 35.95247 44.70905
2023.936 37.84036 34.80260 40.87812 33.19451 42.48621
2023.940 39.32043 36.11723 42.52362 34.42156 44.21929
2023.944 39.28046 35.91994 42.64099 34.14098 44.41994
2023.948 38.45043 34.93963 41.96123 33.08113 43.81974
2023.952 38.70043 35.04553 42.35533 33.11075 44.29012
2023.956 38.92044 35.12691 42.71397 33.11873 44.72215
2023.960 39.54044 35.61316 43.46771 33.53419 45.54668
2023.964 39.72044 35.66383 43.77705 33.51639 45.92448
2023.968 41.09044 36.90849 45.27238 34.69470 47.48617
2023.972 42.59044 38.28680 46.89407 36.00860 49.17228
2023.976 42.84044 38.41846 47.26241 36.07761 49.60326
2023.980 44.15044 39.61321 48.68767 37.21134 51.08953
2023.984 41.92044 37.27081 46.57007 34.80945 49.03143plot_forecast(tecno21_pred,
title = "21 days Close price - Forecast",
Ytitle = "Point Forecast",
Xtitle = "Year")accuracy(tecno21_pred, test28) ME RMSE MAE MPE MAPE MASE
Training set 0.0002181155 0.7664593 0.3852712 -0.01923693 2.450196 0.08150625
Test set 9.9289994958 10.8788967 9.9289995 19.11225159 19.112252 2.10053447
ACF1 Theil's U
Training set 0.001789564 NA
Test set 0.760498080 7.63443Horizonte de 28 días
tecno28_pred <- forecast(auto_model, h = 28)
tecno28_pred Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
2023.904 42.81117 41.78404 43.83830 41.24031 44.38203
2023.908 43.20833 41.75545 44.66122 40.98634 45.43033
2023.912 42.86101 41.12029 44.60174 40.19880 45.52322
2023.916 40.34931 38.33378 42.36484 37.26682 43.43179
2023.920 40.55486 38.29564 42.81408 37.09968 44.01003
2023.924 40.35932 37.88339 42.83524 36.57272 44.14592
2023.928 38.61063 35.93441 41.28686 34.51771 42.70356
2023.932 40.33076 37.46795 43.19357 35.95247 44.70905
2023.936 37.84036 34.80260 40.87812 33.19451 42.48621
2023.940 39.32043 36.11723 42.52362 34.42156 44.21929
2023.944 39.28046 35.91994 42.64099 34.14098 44.41994
2023.948 38.45043 34.93963 41.96123 33.08113 43.81974
2023.952 38.70043 35.04553 42.35533 33.11075 44.29012
2023.956 38.92044 35.12691 42.71397 33.11873 44.72215
2023.960 39.54044 35.61316 43.46771 33.53419 45.54668
2023.964 39.72044 35.66383 43.77705 33.51639 45.92448
2023.968 41.09044 36.90849 45.27238 34.69470 47.48617
2023.972 42.59044 38.28680 46.89407 36.00860 49.17228
2023.976 42.84044 38.41846 47.26241 36.07761 49.60326
2023.980 44.15044 39.61321 48.68767 37.21134 51.08953
2023.984 41.92044 37.27081 46.57007 34.80945 49.03143
2023.988 41.10044 36.34107 45.85981 33.82161 48.37927
2023.992 40.88044 36.01379 45.74708 33.43755 48.32333
2023.996 41.59044 36.61883 46.56204 33.98703 49.19384
2024.000 42.63044 37.55605 47.70482 34.86983 50.39104
2024.004 42.95044 37.77530 48.12557 35.03575 50.86512
2024.008 44.31044 39.03648 49.58439 36.24462 52.37625
2024.012 43.72044 38.34948 49.09139 35.50627 51.93461plot_forecast(tecno28_pred,
title = "28 days Close price - Forecast",
Ytitle = "Point Forecast",
Xtitle = "Year")accuracy(tecno28_pred, test28) ME RMSE MAE MPE MAPE
Training set 2.181155e-04 0.7664593 0.3852712 -0.01923693 2.450196
Test set 1.106789e+01 12.0855244 11.0678896 20.70756616 20.707566
MASE ACF1 Theil's U
Training set 0.08150625 0.001789564 NA
Test set 2.34147294 0.824358743 8.677721Al emplear la función forecast de R para realizar las predicciones, se observa un aumento considerable en la dispersión de las predicciones a medida que aumenta el grado de confianza. Esto es evidente al analizar las métricas de evaluación para diferentes horizontes de predicción.
Para el horizonte de 7 días, la media del error (ME) y la raíz del error cuadrático medio (RMSE) en el conjunto de prueba son significativamente altos, indicando una discrepancia considerable entre las predicciones y los valores observados. Además, los errores porcentuales (MPE y MAPE) son sustancialmente elevados, mostrando una falta de precisión en las predicciones.
A medida que aumenta el horizonte de predicción, esta tendencia de alta dispersión en las predicciones se mantiene. Tanto el RMSE como el MAPE para los horizontes de 14, 21 y 28 días son considerablemente más altos que para el horizonte de 7 días, lo que confirma la tendencia de aumento en la dispersión de las predicciones.
Estos resultados sugieren que la función forecast de R, al no tener en cuenta las medias móviles, puede generar predicciones con una dispersión más alta, especialmente cuando se aumenta el grado de confianza. Para validar la confiabilidad de estas predicciones, es necesario realizar un análisis exhaustivo de los residuos del modelo. Además, se debe considerar la posibilidad de ajustar el modelo para mejorar su desempeño y reducir la dispersión en las predicciones.
Paso 3: Residuales
checkresiduals(fit_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 resultado del test Ljung-Box aplicado a los residuos del modelo ARIMA(3,1,3)(1,1,1)[250] revela hallazgos cruciales sobre la independencia de los residuos, un supuesto fundamental en el análisis de series temporales. Con un valor p sustancialmente menor que el nivel de significancia comúnmente aceptado (p < 2.2e-16), se rechaza la hipótesis nula de independencia de los residuos. Estos resultados implican que los residuos exhiben patrones de correlación significativos en los rezagos, lo que desafía la validez del modelo en términos de independencia. Esta revelación señala la necesidad de reevaluar y posiblemente ajustar el modelo para mejorar la calidad de las predicciones y garantizar la fiabilidad de los resultados obtenidos en el análisis de series temporales.
Paso 4: Criterio BIC y HQIC
Luego de ver los posibles problemas que pueden surgir a partir de la implementación de modelos ARIMA por medio de los anteriores criterios, en esta seccion se procede a estudiar los criterios de BIC y HQIC para encontrar opciones que optimicen las predicciones que se quieran realizar.
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 1Criterio 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 1Al comparar los resultados de los criterios HQIC, BIC y AIC para la selección del modelo óptimo, se observa que los tres criterios coinciden en la especificación del modelo ARIMA(3,1,3)(1,1,1)[250] como el mejor modelo para el conjunto de datos analizado. Esto se refleja en los valores más bajos de los criterios HQIC, BIC y AIC (-7107.385, -7108.267 y -7167.355 respectivamente) asociados con este modelo en comparación con otras especificaciones. Esta consistencia en la selección del modelo sugiere una robustez en la elección del ARIMA(3,1,3)(1,1,1)[250] como el modelo más apropiado para la descripción y predicción de los datos de la serie temporal en cuestión.
##Paso 5: Conclusiones {style=“color: #922B21”}
En esta investigación, se llevó a cabo un exhaustivo análisis de series temporales con el objetivo de modelar y predecir el comportamiento de los datos. Se utilizaron diversas técnicas, incluyendo el método rolling para entrenar modelos con horizontes de predicción de 7, 14, 21 y 28 días, así como la función forecast de R. En la evaluación de los resultados, se observó que las predicciones tendieron a acercarse más a los valores observados en el conjunto de prueba a medida que aumentaba el horizonte de predicción, aunque este patrón varió según el método utilizado.
Al examinar las métricas de evaluación, se notó una mejora gradual en la precisión de las predicciones a medida que se incrementaba el horizonte de predicción. Sin embargo, al utilizar la función forecast de R, que no tiene en cuenta las medias móviles, se observó un aumento significativo en la dispersión de las predicciones a medida que aumentaba el grado de confianza.
Además, se realizó un análisis de los supuestos de normalidad e independencia de los residuos del modelo. Los resultados del test Ljung-Box revelaron que los residuos no eran independientes, lo que cuestiona la validez del modelo en términos de independencia y sugiere la necesidad de realizar ajustes adicionales.
Finalmente, al comparar los criterios de selección de modelos (HQIC, BIC y AIC), se encontró consistencia en la especificación del modelo ARIMA(3,1,3)(1,1,1)[250] como el más apropiado para describir y predecir los datos de la serie temporal analizada.
En conclusión, estos hallazgos resaltan la importancia de una evaluación exhaustiva de los modelos de series temporales, incluyendo la consideración de diferentes métodos de predicción, la validación de supuestos y la selección cuidadosa del modelo final. Este proceso garantiza la obtención de predicciones confiables y útiles para la toma de decisiones en diversos contextos.