Considere la serie de tiempo asociada con las acciones de Tecnoglass desde que comenzó a comercializarse hasta la fecha del día de hoy. Puede utilizar la API de Yahoo Finance para obtener esta serie de tiempo (ver yahoofinancer).
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. Cualquier cantidad de pasos extra para enriquecer su análisis predictivo serán aceptados siempre y cuando sean acordes con lo que indica la teoría de análisis de series de tiempo.
Repita el paso 2 ahora sin utilizar rolling. Esto es, realice el pronóstico solo utilizando forecast() para los diferentes horizontes de predicción, 7, 14 días, 21 días, 28 días.
Realice tablas de error para los ítems 1 y 2, utilizando las métricas: MAPE, MAE, RMSE, MSE, R2. Además, agregue el gráfico de correlación entre la observación real y su predicción en el test, Corr(yt,yt)Corr(��,��).
Repita el análisis desarrollado en los pasos anteriores, considerando ahora el criterio de inferencia Bayesiana (BIC) y el criterio de información de Hannan–Quinn (HQIC) para encontrar el mejor modelo ARIMA y, compare los errores con aquellos obtenidos con el criterio de Akaike.
Escriba en cada paso las conclusiones y análisis estadísticos asociados con los resultados obtenidos. Realice tests de normalidad e independencia para los residuales obtenidos para cada predicción, en cada caso agregue las correspondientes conclusiones. Figuras y algoritmos que no estén acompañados de una conclusión, descripción y análisis estadístico, no serán tenidas en cuenta.
MODELOS ARIMA - Minimizando AIC
Vamos a encontrar ahora el mejor orden de los parámetros p,q,d con una función que minimiza el AIC
best_aic <- Inf
best_pdq <- NULL
best_PDQ <- NULL
fit <- NULL
p_n<-3
d_n<-3
q_n<-2
for(p in 1:p_n) {
print(paste("Iniciando en p =", p))
for(d in 1:d_n) {
print(paste("iniciando en d= ", d))
for (q in 1:q_n) {
print(paste("iniciando en q=", q))
for(P in 1:p_n) {
for(D in 1:d_n) {
for (Q in 1:q_n) {
tryCatch({
fit <- arima(scale(ts),
order=c(p, d, q),
seasonal = list(order = c(P, D, Q), period = 12),
xreg=1:length(ts),
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)
}
print(best_aic)
print(best_pdq)
}, error=function(e){})
}
}
}
}
}
}
Ahora reentrenamos el modelo usando los parámetros encontrados en el paso anterior:
Capturamos y evaluamos los residuales:
Funcion para evaluar residuales
evaluar_residuales <- function(model, alpha = 0.05) {
# Prueba de normalidad (Shapiro-Wilk)
shapiro_test <- shapiro.test(model$residuals)
if (shapiro_test$p.value > alpha) {
resultado_normalidad <- "Los residuales se distribuyen normalmente"
} else {
resultado_normalidad <- "Los residuales no se distribuyen normalmente"
}
# Prueba de independencia (Ljung-Box)
ljung <- Box.test(model$residuals, type = "Ljung-Box", lag = 20)
ljung_pvalues <- ljung$p.value
if (all(ljung_pvalues > alpha)) {
resultado_independencia <- "Los residuales son independientes"
} else {
resultado_independencia <- "Los residuales no son independientes"
}
# Gráfico de residuos
checkresiduals(model)
return(list(resultado_normalidad = resultado_normalidad, resultado_independencia = resultado_independencia))
}
best_fit <- arima(scale(ts), order=c(2,1,2), seasonal = list(order = c(1,1,1), period = 12), xreg=1:length(ts), method="CSS-ML")
best_residuales<- residuals(best_fit)
evaluar_residuales(best_fit,alpha = 0.05)
Ljung-Box test
data: Residuals from ARIMA(2,1,2)(1,1,1)[12]
Q* = 1054.7, df = 724, p-value = 9.77e-15
Model df: 6. Total lags used: 730
$resultado_normalidad
[1] "Los residuales no se distribuyen normalmente"
$resultado_independencia
[1] "Los residuales son independientes"
Aunque los residuales no se distribuyan normalmente, podemos confirmar que son independientes lo cual nos garantiza que el modelo esta capturando correctamente la correlación de la serie y el los residuales son ruido
Funciones
Arima Rolling
arima_rolling <- function(history, test, best_order) {
predictions <- vector()
residuals <- vector()
for (t in 1:length(test)) {
model <- arima(history, order = best_order)
model_fit <- forecast:::Arima(y = history, model = model, h = 1)
residuals <- residuals + model$residuals
output <- forecast:::forecast.Arima(model_fit, h = 1)
yhat <- output$mean[1]
predictions <- c(predictions, yhat)
obs <- test[t]
history <- c(history, obs)
cat(sprintf("predicted=%f, expected=%f\n", yhat, obs))
}
return(list(predictions = predictions, residuals = residuals))
}
Arima sin Rolling
arima_sin_rolling <- function(train, test, best_order) {
model <- arima(train, order = best_order) # Aquí se utiliza 'train' en lugar de 'history'
predictions <- forecast::forecast(model, h = length(test))$mean
residuals <- residuals(model)
return(list(predictions = predictions, model=model))
}
Graficos de predicciones
# Función para graficar predicciones de forma interactiva con plotly
grafico_predicciones<- function(dates_train, train, dates_w, test_wl, yhat_w, window) {
# Crear el gráfico de líneas con plotly
p <- plot_ly() %>%
add_lines(x = dates_train, y = train, name = "Train", color = I("green")) %>%
add_lines(x = dates_w, y = test_wl, name = "Test", color = I("blue")) %>%
add_lines(x = dates_w, y = yhat_w, name = "Forecast", color = I("red")) %>%
layout(title = paste("Gráfico de Predicciones para Horizonte de", window, "Días"),
xaxis = list(title = "Fecha"),
yaxis = list(title = "Valor"),
legend = list(x = 0.9, y = 1)
)
# Mostrar el gráfico
p
}
Metricas del modelo
# Función para calcular las métricas
metrics <- function(forecast, actual, str_name) {
mape <- mean(abs((forecast - actual) / actual)) # MAPE
mae <- mean(abs(forecast - actual)) # MAE
rmse <- sqrt(mean((forecast - actual)^2)) # RMSE
mse <- mean((forecast - actual)^2) # MSE
r2 <- cor(forecast, actual)^2 # R^2
df_metrics <- data.frame(MAE = mae,
MSE = mse,
MAPE = mape,
RMSE = rmse,
R2 = r2,
row.names = str_name)
return(kable(df_metrics, caption = "Métricas de Predicción"))
}
Predicciones usando ARIMAS, métricas y gráficos
Para ventana de 7 días
n_BTC <- length(ts)
n_test <- 7
train_size <- (n_BTC - n_test)
cat("1. FRACCIONAMOS DATASET--------\n")
1. FRACCIONAMOS DATASET--------
# Se asume que la fecha de inicio ya está definida como start_date
print(paste("La fecha inicial es: ", start_date))
[1] "La fecha inicial es: 2014-04-01"
end_date <- start_date + (train_size - 1)
print(paste("La fecha final es: ", end_date))
[1] "La fecha final es: 2024-03-22"
train <- ts[1:train_size]
print(paste("Train es: ",length (train)))
test_w <- ts[train_size:(train_size + n_test-1)]
print(paste("Test es: ", length(test_w)))
# Creamos la secuencia de fechas de entrenamiento y prueba
dates_train <- seq(start_date, end_date, by = "days")
print(paste("Fechas de entrenamiento:", length(dates_train)))
[1] "Fechas de entrenamiento: 3644"
start_date_test <- end_date + 1
end_date_test <- start_date_test + n_test-1
dates_test <- seq(start_date_test, end_date_test, by = "days")
print(paste("Fechas de prueba:", length(dates_test)))
[1] "Fechas de prueba: 7"
ARIMA CON ROLLING
predictions <- vector()
residuals <- vector()
model <- arima(train, order = best_pdq)
history<-train
for (t in 1:length(test_w)) {
# Realizar el pronóstico para un paso adelante
model_fit <- forecast:::Arima(y = history, model = model, h = 1)
# Actualizar los residuales
residuals <- c(residuals, model_fit$residuals)
# Obtener la predicción
output <- forecast:::forecast.Arima(model_fit, h = 1)
yhat <- output$mean[1]
# Almacenar la predicción
predictions <- c(predictions, yhat)
# Avanzar al siguiente paso de tiempo
obs <- test_w[t]
history <- c(history, obs) # Agregar el nuevo valor observado a history
history <- history[-1] # Eliminar el primer valor de history para mantener su longitud constante
# Imprimir la predicción y el valor esperado
cat(sprintf("predicted=%f, expected=%f\n", yhat, obs))
}
predicted=9.524254, expected=9.520000
predicted=9.517566, expected=9.500000
predicted=9.502520, expected=9.490000
predicted=9.487774, expected=9.420000
predicted=9.422357, expected=9.500000
predicted=9.498389, expected=9.450000
predicted=9.450816, expected=9.480000
# Devolver las predicciones y los residuales en una lista
evaluar_residuales(model_fit,alpha = 0.05)
Ljung-Box test
data: Residuals from ARIMA(2,1,2)
Q* = 11.121, df = 6, p-value = 0.08471
Model df: 4. Total lags used: 10
$resultado_normalidad
[1] "Los residuales no se distribuyen normalmente"
$resultado_independencia
[1] "Los residuales son independientes"
grafico_predicciones(dates_train,train,dates_test, test_w, predictions, n_test)
metrics(predictions,test_w,"7_days")
Métricas de Predicción
| 7_days |
0.0367613 |
0.0020426 |
0.0038831 |
0.0451953 |
0.0057914 |
list(predictions = predictions)
$predictions
[1] 9.524254 9.517566 9.502520 9.487774 9.422357 9.498389 9.450816
ARIMA SIN ROLLING
resultados <- arima_sin_rolling(train, test_w, best_pdq)
predictions <- as.numeric(resultados$predictions)
modelo <- resultados$model
evaluar_residuales(modelo, alpha = 0.05)
Ljung-Box test
data: Residuals from ARIMA(2,1,2)
Q* = 11.119, df = 6, p-value = 0.08476
Model df: 4. Total lags used: 10
$resultado_normalidad
[1] "Los residuales no se distribuyen normalmente"
$resultado_independencia
[1] "Los residuales son independientes"
metrics(predictions,test_w,"7_days")
Métricas de Predicción
| 7_days |
0.0432231 |
0.002824 |
0.0045704 |
0.0531415 |
0.3947983 |
Para ventana de 14 días
ARIMA CON ROLLING
n_BTC <- length(ts)
n_test <- 14
train_size <- (n_BTC - n_test)
cat("1. FRACCIONAMOS DATASET--------\n")
1. FRACCIONAMOS DATASET--------
# Se asume que la fecha de inicio ya está definida como start_date
print(paste("La fecha inicial es: ", start_date))
[1] "La fecha inicial es: 2014-04-01"
end_date <- start_date + (train_size - 1)
print(paste("La fecha final es: ", end_date))
[1] "La fecha final es: 2024-03-15"
train <- ts[1:train_size]
print(paste("Train es: ",length (train)))
test_w <- ts[train_size:(train_size + n_test-1)]
print(paste("Test es: ", length(test_w)))
# Creamos la secuencia de fechas de entrenamiento y prueba
dates_train <- seq(start_date, end_date, by = "days")
print(paste("Fechas de entrenamiento:", length(dates_train)))
[1] "Fechas de entrenamiento: 3637"
start_date_test <- end_date + 1
end_date_test <- start_date_test + n_test-1
dates_test <- seq(start_date_test, end_date_test, by = "days")
print(paste("Fechas de prueba:", length(dates_test)))
[1] "Fechas de prueba: 14"
predictions <- vector()
residuals <- vector()
model <- arima(train, order = best_pdq)
history<-train
for (t in 1:length(test_w)) {
# Realizar el pronóstico para un paso adelante
model_fit <- forecast:::Arima(y = history, model = model, h = 1)
# Actualizar los residuales
residuals <- c(residuals, model_fit$residuals)
# Obtener la predicción
output <- forecast:::forecast.Arima(model_fit, h = 1)
yhat <- output$mean[1]
# Almacenar la predicción
predictions <- c(predictions, yhat)
# Avanzar al siguiente paso de tiempo
obs <- test_w[t]
history <- c(history, obs) # Agregar el nuevo valor observado a history
history <- history[-1] # Eliminar el primer valor de history para mantener su longitud constante
# Imprimir la predicción y el valor esperado
cat(sprintf("predicted=%f, expected=%f\n", yhat, obs))
}
predicted=8.954019, expected=8.950000
predicted=8.945793, expected=9.040000
predicted=9.044093, expected=9.150000
predicted=9.145080, expected=9.740000
predicted=9.743877, expected=9.420000
predicted=9.408154, expected=9.570000
predicted=9.581364, expected=9.690000
predicted=9.677175, expected=9.520000
predicted=9.529908, expected=9.500000
predicted=9.490955, expected=9.490000
predicted=9.499388, expected=9.420000
predicted=9.411145, expected=9.500000
predicted=9.509772, expected=9.450000
predicted=9.440020, expected=9.480000
# Devolver las predicciones y los residuales en una lista
evaluar_residuales(model_fit,alpha = 0.05)
Ljung-Box test
data: Residuals from ARIMA(2,1,2)
Q* = 11.07, df = 6, p-value = 0.08625
Model df: 4. Total lags used: 10
$resultado_normalidad
[1] "Los residuales no se distribuyen normalmente"
$resultado_independencia
[1] "Los residuales son independientes"
metrics(predictions,test_w,"14_days")
Métricas de Predicción
| 14_days |
0.1321032 |
0.0401355 |
0.0138651 |
0.2003385 |
0.4276038 |
list(predictions = predictions)
$predictions
[1] 8.954019 8.945793 9.044093 9.145080 9.743877 9.408154 9.581364 9.677175
[9] 9.529908 9.490955 9.499388 9.411145 9.509772 9.440020
ARIMA SIN ROLLING
resultados <- arima_sin_rolling(train, test_w, best_pdq)
predictions <- resultados$predictions
modelo <- resultados$model
evaluar_residuales(modelo, alpha = 0.05)
Ljung-Box test
data: Residuals from ARIMA(2,1,2)
Q* = 11.074, df = 6, p-value = 0.08611
Model df: 4. Total lags used: 10
$resultado_normalidad
[1] "Los residuales no se distribuyen normalmente"
$resultado_independencia
[1] "Los residuales son independientes"
metrics(predictions,test_w,"14_days")
Métricas de Predicción
| 14_days |
0.4716584 |
0.2698728 |
0.049535 |
0.5194928 |
0.0930458 |
Para la ventana de 21 días
ARIMA CON ROLLING
n_BTC <- length(ts)
n_test <- 21
train_size <- (n_BTC - n_test)
cat("1. FRACCIONAMOS DATASET--------\n")
1. FRACCIONAMOS DATASET--------
# Se asume que la fecha de inicio ya está definida como start_date
print(paste("La fecha inicial es: ", start_date))
[1] "La fecha inicial es: 2014-04-01"
end_date <- start_date + (train_size - 1)
print(paste("La fecha final es: ", end_date))
[1] "La fecha final es: 2024-03-08"
train <- ts[1:train_size]
print(paste("Train es: ",length (train)))
test_w <- ts[train_size:(train_size + n_test-1)]
print(paste("Test es: ", length(test_w)))
# Creamos la secuencia de fechas de entrenamiento y prueba
dates_train <- seq(start_date, end_date, by = "days")
print(paste("Fechas de entrenamiento:", length(dates_train)))
[1] "Fechas de entrenamiento: 3630"
start_date_test <- end_date + 1
end_date_test <- start_date_test + n_test-1
dates_test <- seq(start_date_test, end_date_test, by = "days")
print(paste("Fechas de prueba:", length(dates_test)))
[1] "Fechas de prueba: 21"
predictions <- vector()
residuals <- vector()
model <- arima(train, order = best_pdq)
history<-train
for (t in 1:length(test_w)) {
# Realizar el pronóstico para un paso adelante
model_fit <- forecast:::Arima(y = history, model = model, h = 1)
# Actualizar los residuales
residuals <- c(residuals, model_fit$residuals)
# Obtener la predicción
output <- forecast:::forecast.Arima(model_fit, h = 1)
yhat <- output$mean[1]
# Almacenar la predicción
predictions <- c(predictions, yhat)
# Avanzar al siguiente paso de tiempo
obs <- test_w[t]
history <- c(history, obs) # Agregar el nuevo valor observado a history
history <- history[-1] # Eliminar el primer valor de history para mantener su longitud constante
# Imprimir la predicción y el valor esperado
cat(sprintf("predicted=%f, expected=%f\n", yhat, obs))
}
predicted=8.810572, expected=8.820000
predicted=8.828913, expected=8.690000
predicted=8.680861, expected=8.980000
predicted=8.990808, expected=9.050000
predicted=9.037239, expected=8.950000
predicted=8.960137, expected=8.850000
predicted=8.839737, expected=8.950000
predicted=8.961201, expected=8.950000
predicted=8.938561, expected=9.040000
predicted=9.051060, expected=9.150000
predicted=9.138294, expected=9.740000
predicted=9.750801, expected=9.420000
predicted=9.401416, expected=9.570000
predicted=9.587580, expected=9.690000
predicted=9.670969, expected=9.520000
predicted=9.535719, expected=9.500000
predicted=9.484969, expected=9.490000
predicted=9.505076, expected=9.420000
predicted=9.405490, expected=9.500000
predicted=9.515285, expected=9.450000
predicted=9.434673, expected=9.480000
# Devolver las predicciones y los residuales en una lista
evaluar_residuales(model_fit,alpha = 0.05)
Ljung-Box test
data: Residuals from ARIMA(2,1,2)
Q* = 11.091, df = 6, p-value = 0.0856
Model df: 4. Total lags used: 10
$resultado_normalidad
[1] "Los residuales no se distribuyen normalmente"
$resultado_independencia
[1] "Los residuales son independientes"
metrics(predictions,test_w,"21_days")
Métricas de Predicción
| 21_days |
0.1291102 |
0.0343545 |
0.0138381 |
0.1853498 |
0.6953025 |
list(predictions = predictions)
$predictions
[1] 8.810572 8.828913 8.680861 8.990808 9.037239 8.960137 8.839737 8.961201
[9] 8.938561 9.051060 9.138294 9.750801 9.401416 9.587580 9.670969 9.535719
[17] 9.484969 9.505076 9.405490 9.515285 9.434673