Series de Tiempo
Maria Paula Alvarez y Daniel Piñeres
2024-04-14
Tecnoglass
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).
library(quantmod)## Warning: package 'quantmod' was built under R version 4.3.3## Loading required package: xts## Warning: package 'xts' was built under R version 4.3.3## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: TTR## Warning: package 'TTR' was built under R version 4.3.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(plotly)## Warning: package 'plotly' was built under R version 4.3.3## Loading required package: ggplot2##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutlibrary(ggplot2)
library(forecast)## Warning: package 'forecast' was built under R version 4.3.3# Símbolo de Tecnoglass en Yahoo Finance
ticker <- "TGLS"
# Descarga los datos históricos desde una fecha específica hasta hoy
start_date <- as.Date("2016-04-13")
end_date <- Sys.Date() -1
# Obtén los datos
getSymbols(ticker, src = "yahoo", from = start_date, to = end_date)## [1] "TGLS"df_TGLS= data.frame(TGLS)
df_TGLS$Date <- rownames(df_TGLS)
df_TGLS$Date <- as.Date(df_TGLS$Date, format = "%Y-%m-%d")
plot_ly(data = df_TGLS, x = ~Date, y = ~TGLS.Close, type = "scatter", mode = "lines",
text = ~paste("Fecha: ", Date, "<br>Precio de cierre: $", TGLS.Close)) %>%
layout(title = "Gráfico interactivo de Precio de Cierre de TGLS",
xaxis = list(title = "Fecha"),
yaxis = list(title = "Precio de cierre"))Se observa que la serie de tiempo no es estacionaria, al evidenciarse tendencia y varianza no constante.A partir de mayo del 2021 se evidencia un alza en los precios de cierre de las acciones de Tecnoglass, el COO de Tecnoglass anunció que abrirían una nueva planta en mayo con una inversión de US$40 millones. Esta planta se consideraría la fábrica de vidrio “más grande y moderna del mundo”.
## Warning: package 'TSstudio' was built under R version 4.3.3# Check the information of the new ts object
ts_info(TGLS)## The TGLS series is a xts object with 6 variables and 2015 observations
## Frequency: daily
## Start time: 2016-04-13
## End time: 2024-04-15library(xts)
TGLS_xts <- xts(TGLS$TGLS.Close, frequency = 365)
# Check the information of the new ts object
ts_info(TGLS_xts)## The TGLS_xts series is a xts object with 1 variable and 2015 observations
## Frequency: daily
## Start time: 2016-04-13
## End time: 2024-04-15# Check the information of the new ts object
TGLS_ts = ts(TGLS_xts, frequency = 365)
ts_decompose(TGLS_ts)En la gráfica se observa un patrón estacional adicional a la tendencia identificada previamente.
library(tseries)## Warning: package 'tseries' was built under R version 4.3.3adf_result <- adf.test(TGLS_ts)
# Imprimir los resultados
cat("ADF Statistic:", adf_result$statistic, "\n")## ADF Statistic: -0.9003795cat("p-value:", adf_result$p.value, "\n")## p-value: 0.9526583# Interpretación de los resultados
if (adf_result$p.value < 0.05) {
cat("Conclusión: Rechazamos la hipótesis nula. La serie temporal es estacionaria.\n")
} else {
cat("Conclusión: No rechazamos la hipótesis nula. La serie temporal es no estacionaria.\n")
}## Conclusión: No rechazamos la hipótesis nula. La serie temporal es no estacionaria.TGLS_split <- ts_split(ts(TGLS$TGLS.Close), sample.out = 28)
train <- TGLS_split$train
test <- TGLS_split$testSe realiza la separación de los datos para entrenamiento y testing.
par(mfrow=c(1,2))
acf(train, lag.max = 60)
pacf(train, lag.max = 60)
Como la serie es no estacionaria, procederemos con la
diferenciación.
En la gráfica proporcionada, el número de retrasos significativos en el PACF (Función de Autocorrelación Parcial) se puede determinar observando los picos que exceden el umbral de significancia. En este caso, parece que el primer pico significativo está en el retraso 1.
La autocorrelación de la serie con sus rezagos decae lentamente en el tiempo de forma lineal. La eliminación de la tendencia de la serie y de la correlación entre la serie y sus rezagos puede hacerse por diferenciación.
par(mfrow=c(1,3))
acf(diff(TGLS_ts, lag = 1), lag.max = 60)
pacf(diff(TGLS_ts, lag = 1), lag.max = 60)
ts_plot(diff(TGLS_xts, lag = 1),
title = "First Differencing",
Xtitle = "Year",
Ytitle = "Differencing Close Tecnoglass")
adf_result <- adf.test(diff(TGLS_ts, lag = 1))## Warning in adf.test(diff(TGLS_ts, lag = 1)): p-value smaller than printed
## p-value# Imprimir los resultados
cat("ADF Statistic:", adf_result$statistic, "\n")## ADF Statistic: -12.35077cat("p-value:", adf_result$p.value, "\n")## p-value: 0.01# Interpretación de los resultados
if (adf_result$p.value < 0.05) {
cat("Conclusión: Rechazamos la hipótesis nula. La serie temporal es estacionaria.\n")
} else {
cat("Conclusión: No rechazamos la hipótesis nula. La serie temporal es no estacionaria.\n")
}## Conclusión: Rechazamos la hipótesis nula. La serie temporal es estacionaria.2. Repita TODOS los pasos indicados en esta sección para encontrar modelos ARIMA para predecir el precio de Bitcoin con los siguientes horizontes: 7, 14, 21 y 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.
best_aic <- Inf
best_order <- NULL
best_mdl <- NULL
pq_rng <- seq(0, 7)
d_rng <- seq(0, 2)
for (i in pq_rng) {
for (d in d_rng) {
for (j in pq_rng) {
tryCatch({
tmp_mdl <- arima(train, order = c(i, d, j))
tmp_aic <- AIC(tmp_mdl)
if (tmp_aic < best_aic) {
best_aic <- tmp_aic
best_order <- c(i, d, j)
best_mdl <- tmp_mdl
}
}, error = function(e) {})
}
}
}## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1## Warning in log(s2): NaNs produced## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1## Warning in log(s2): NaNs produced
## Warning in log(s2): NaNs produced
## Warning in log(s2): NaNs produced## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1## Warning in log(s2): NaNs produced## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1print(best_order)## [1] 5 1 7best_order <- c(5,1,7)El mejor modelo ARIMA basado en el AIC un ARIMA utiliza los últimos 5 valores anteriores, aplica una diferenciación de grado 1 y utiliza una media móvil basada en los últimos 7 datos para hacer predicciones sobre la serie temporal.
Rolling forecasting para predecir cifras futuras de forma continua durante un periodo de tiempo, este caso se utilizaran los periodos 7, 14,21 y 28 días.
forecast_accuracy <- function(forecast, actual) {
mape <- mean(abs((forecast - actual) / actual))
mae <- mean(abs(forecast - actual))
rmse <- sqrt(mean((forecast - actual)^2))
mse <- mean((forecast - actual)^2)
r2 <- 1 - sum((actual - forecast)^2) / sum((actual - mean(actual))^2)
return(list(MAPE = mape, MAE = mae, RMSE = rmse, MSE = mse, R2 = r2))
}non_rolling_arima_forecast <- function(data, horizons, order) {
results <- list()
residuals <- list()
for (horizon in horizons) {
train_size <- length(data) - horizon
train <- data[1:train_size]
test <- data[(train_size + 1):length(data)]
model <- Arima(train, order = order)
forecast_result <- forecast(model, h = horizon)
# Calculate residuals and forecast accuracy
forecast_length <- length(forecast_result$mean)
if (forecast_length == length(test)) {
current_residuals <- forecast_result$mean - test
residuals[[as.character(horizon)]] <- current_residuals
accuracy <- forecast_accuracy(forecast_result$mean, test)
results[[as.character(horizon)]] <- list(Forecast = forecast_result$mean, Accuracy = accuracy)
} else {
stop("Forecast and test data lengths do not match.")
}
}
return(list(Results = results, Residuals = residuals))
}TGLS_ts = ts(TGLS_xts, frequency = 1)
# Define los horizontes
horizons <- c(7, 14, 21, 28)
# Ejecutar el modelo ARIMA para cada horizonte y evaluar
non_rolling_results <- non_rolling_arima_forecast(TGLS_ts, horizons, best_order)non_rolling_results## $Results
## $Results$`7`
## $Results$`7`$Forecast
## Time Series:
## Start = 2009
## End = 2015
## Frequency = 1
## [1] 52.50531 52.54140 52.61208 52.80297 52.60863 52.30420 52.20950
##
## $Results$`7`$Accuracy
## $Results$`7`$Accuracy$MAPE
## [1] 0.0850605
##
## $Results$`7`$Accuracy$MAE
## [1] 4.942273
##
## $Results$`7`$Accuracy$RMSE
## [1] 5.271315
##
## $Results$`7`$Accuracy$MSE
## [1] 27.78676
##
## $Results$`7`$Accuracy$R2
## [1] -7.190476
##
##
##
## $Results$`14`
## $Results$`14`$Forecast
## Time Series:
## Start = 2002
## End = 2015
## Frequency = 1
## [1] 53.49365 52.98504 52.95737 53.05907 52.81746 52.57484 52.92077 53.23539
## [9] 53.14952 53.11159 53.41041 53.34143 52.91586 52.76168
##
## $Results$`14`$Accuracy
## $Results$`14`$Accuracy$MAPE
## [1] 0.04458097
##
## $Results$`14`$Accuracy$MAE
## [1] 2.545381
##
## $Results$`14`$Accuracy$RMSE
## [1] 3.411227
##
## $Results$`14`$Accuracy$MSE
## [1] 11.63647
##
## $Results$`14`$Accuracy$R2
## [1] -0.3136622
##
##
##
## $Results$`21`
## $Results$`21`$Forecast
## Time Series:
## Start = 1995
## End = 2015
## Frequency = 1
## [1] 45.04813 44.83100 45.08227 45.15161 44.97562 45.11003 45.50195 45.45809
## [9] 45.19023 45.22012 45.31910 45.02171 44.81268 45.03856 45.25340 45.13310
## [17] 45.19183 45.47889 45.47455 45.14484 45.08242
##
## $Results$`21`$Accuracy
## $Results$`21`$Accuracy$MAPE
## [1] 0.1452658
##
## $Results$`21`$Accuracy$MAE
## [1] 7.959041
##
## $Results$`21`$Accuracy$RMSE
## [1] 8.818567
##
## $Results$`21`$Accuracy$MSE
## [1] 77.76713
##
## $Results$`21`$Accuracy$R2
## [1] -4.140183
##
##
##
## $Results$`28`
## $Results$`28`$Forecast
## Time Series:
## Start = 1988
## End = 2015
## Frequency = 1
## [1] 44.16424 44.28936 44.69218 44.61214 44.41224 44.62790 44.64346 44.26862
## [9] 44.08936 44.32586 44.33465 44.21355 44.40112 44.72411 44.62934 44.45628
## [17] 44.50077 44.49294 44.19811 44.11319 44.30078 44.40809 44.35543 44.52022
## [25] 44.70024 44.61510 44.40358 44.39545
##
## $Results$`28`$Accuracy
## $Results$`28`$Accuracy$MAPE
## [1] 0.1227172
##
## $Results$`28`$Accuracy$MAE
## [1] 6.674346
##
## $Results$`28`$Accuracy$RMSE
## [1] 8.244068
##
## $Results$`28`$Accuracy$MSE
## [1] 67.96465
##
## $Results$`28`$Accuracy$R2
## [1] -1.865966
##
##
##
##
## $Residuals
## $Residuals$`7`
## Time Series:
## Start = 2009
## End = 2015
## Frequency = 1
## [1] -1.434692 -2.898600 -5.797922 -6.197033 -6.391369 -5.365794 -6.510503
##
## $Residuals$`14`
## Time Series:
## Start = 2002
## End = 2015
## Frequency = 1
## [1] -0.106352126 0.305039298 0.927369637 1.509073840 2.067458783
## [6] -0.005164509 0.460767218 -0.704611842 -2.290474001 -5.298409707
## [11] -5.589589827 -5.658565941 -4.754134071 -5.958325780
##
## $Residuals$`21`
## Time Series:
## Start = 1995
## End = 2015
## Frequency = 1
## [1] -0.7518699 -0.5889992 -2.4677281 -5.2783915 -7.9243860 -7.1499664
## [7] -7.9680475 -8.1419073 -7.4897705 -6.8098741 -6.2308976 -5.7282908
## [13] -7.7673199 -7.4214364 -8.6865948 -10.3069033 -13.2181732 -13.5211087
## [19] -13.5254472 -12.5251587 -13.6375818
##
## $Residuals$`28`
## Time Series:
## Start = 1988
## End = 2015
## Frequency = 1
## [1] -0.36575852 -0.72063844 -0.18782047 -0.03786115 -0.58775641
## [6] -0.95209855 -0.81654404 -1.53137743 -1.33064031 -3.22413934
## [11] -6.09535248 -8.68645419 -7.85887488 -8.74589012 -8.97065909
## [16] -8.22372471 -7.52923335 -7.05706307 -6.55188689 -8.46681427
## [21] -8.15921676 -9.53190402 -11.08456744 -13.88977846 -14.29975882
## [26] -14.38489809 -13.26641938 -14.32455094rolling_arima_forecast <- function(data, horizons, order) {
results <- list()
residuals <- list()
for (horizon in horizons) {
train_size <- length(data) - horizon
train <- window(data, end = train_size)
test <- window(data, start = train_size + 1)
predictions <- vector(mode = "numeric", length = length(test))
for (i in 1:length(test)) {
model <- Arima(train, order = order)
forecast_result <- forecast(model, h = 1)
prediction <- as.numeric(forecast_result$mean)
# Assign predictions to the appropriate index
predictions[i] <- prediction
}
predictions_ts <- ts(predictions, start = start(test), frequency = frequency(test))
# Calculate residuals and forecast accuracy
if (length(predictions_ts) == length(test)) {
current_residuals <- predictions_ts - test
residuals[[as.character(horizon)]] <- current_residuals
accuracy <- forecast_accuracy(predictions_ts, test)
results[[as.character(horizon)]] <- list(Forecast = predictions_ts, Accuracy = accuracy)
} else {
stop("Predictions and test data lengths do not match.")
}
}
return(list(Results = results, Residuals = residuals))
}TGLS_ts = ts(TGLS_xts, frequency = 1)
# Define los horizontes
horizons <- c(7, 14, 21, 28)
# Ejecutar el modelo ARIMA para cada horizonte y evaluar
rolling_results <- rolling_arima_forecast(TGLS_ts, horizons, best_order)rolling_results## $Results
## $Results$`7`
## $Results$`7`$Forecast
## Time Series:
## Start = 2009
## End = 2015
## Frequency = 1
## [1] 52.50531 52.50531 52.50531 52.50531 52.50531 52.50531 52.50531
##
## $Results$`7`$Accuracy
## $Results$`7`$Accuracy$MAPE
## [1] 0.08516772
##
## $Results$`7`$Accuracy$MAE
## [1] 4.948978
##
## $Results$`7`$Accuracy$RMSE
## [1] 5.280621
##
## $Results$`7`$Accuracy$MSE
## [1] 27.88496
##
## $Results$`7`$Accuracy$R2
## [1] -7.21942
##
##
##
## $Results$`14`
## $Results$`14`$Forecast
## Time Series:
## Start = 2002
## End = 2015
## Frequency = 1
## [1] 53.49365 53.49365 53.49365 53.49365 53.49365 53.49365 53.49365 53.49365
## [9] 53.49365 53.49365 53.49365 53.49365 53.49365 53.49365
##
## $Results$`14`$Accuracy
## $Results$`14`$Accuracy$MAPE
## [1] 0.04643121
##
## $Results$`14`$Accuracy$MAE
## [1] 2.624479
##
## $Results$`14`$Accuracy$RMSE
## [1] 3.268669
##
## $Results$`14`$Accuracy$MSE
## [1] 10.6842
##
## $Results$`14`$Accuracy$R2
## [1] -0.206158
##
##
##
## $Results$`21`
## $Results$`21`$Forecast
## Time Series:
## Start = 1995
## End = 2015
## Frequency = 1
## [1] 45.04813 45.04813 45.04813 45.04813 45.04813 45.04813 45.04813 45.04813
## [9] 45.04813 45.04813 45.04813 45.04813 45.04813 45.04813 45.04813 45.04813
## [17] 45.04813 45.04813 45.04813 45.04813 45.04813
##
## $Results$`21`$Accuracy
## $Results$`21`$Accuracy$MAPE
## [1] 0.1473899
##
## $Results$`21`$Accuracy$MAE
## [1] 8.078537
##
## $Results$`21`$Accuracy$RMSE
## [1] 8.966159
##
## $Results$`21`$Accuracy$MSE
## [1] 80.39201
##
## $Results$`21`$Accuracy$R2
## [1] -4.31368
##
##
##
## $Results$`28`
## $Results$`28`$Forecast
## Time Series:
## Start = 1988
## End = 2015
## Frequency = 1
## [1] 44.16424 44.16424 44.16424 44.16424 44.16424 44.16424 44.16424 44.16424
## [9] 44.16424 44.16424 44.16424 44.16424 44.16424 44.16424 44.16424 44.16424
## [17] 44.16424 44.16424 44.16424 44.16424 44.16424 44.16424 44.16424 44.16424
## [25] 44.16424 44.16424 44.16424 44.16424
##
## $Results$`28`$Accuracy
## $Results$`28`$Accuracy$MAPE
## [1] 0.1278024
##
## $Results$`28`$Accuracy$MAE
## [1] 6.934688
##
## $Results$`28`$Accuracy$RMSE
## [1] 8.473741
##
## $Results$`28`$Accuracy$MSE
## [1] 71.80429
##
## $Results$`28`$Accuracy$R2
## [1] -2.027878
##
##
##
##
## $Residuals
## $Residuals$`7`
## Time Series:
## Start = 2009
## End = 2015
## Frequency = 1
## test
## [1,] -1.434692
## [2,] -2.934692
## [3,] -5.904693
## [4,] -6.494693
## [5,] -6.494693
## [6,] -5.164691
## [7,] -6.214694
##
## $Residuals$`14`
## Time Series:
## Start = 2002
## End = 2015
## Frequency = 1
## test
## [1,] -0.1063521
## [2,] 0.8136460
## [3,] 1.4636476
## [4,] 1.9436471
## [5,] 2.7436463
## [6,] 0.9136445
## [7,] 1.0336473
## [8,] -0.4463523
## [9,] -1.9463523
## [10,] -4.9163535
## [11,] -5.5063537
## [12,] -5.5063537
## [13,] -4.1763518
## [14,] -5.2263549
##
## $Residuals$`21`
## Time Series:
## Start = 1995
## End = 2015
## Frequency = 1
## test
## [1,] -0.7518699
## [2,] -0.3718689
## [3,] -2.5018699
## [4,] -5.3818710
## [5,] -7.8518722
## [6,] -7.2118690
## [7,] -8.4218719
## [8,] -8.5518692
## [9,] -7.6318710
## [10,] -6.9818695
## [11,] -6.5018699
## [12,] -5.7018707
## [13,] -7.5318725
## [14,] -7.4118698
## [15,] -8.8918693
## [16,] -10.3918693
## [17,] -13.3618705
## [18,] -13.9518707
## [19,] -13.9518707
## [20,] -12.6218689
## [21,] -13.6718719
##
## $Residuals$`28`
## Time Series:
## Start = 1988
## End = 2015
## Frequency = 1
## test
## [1,] -0.3657585
## [2,] -0.8457581
## [3,] -0.7157608
## [4,] -0.4857613
## [5,] -0.8357597
## [6,] -1.4157616
## [7,] -1.2957588
## [8,] -1.6357590
## [9,] -1.2557579
## [10,] -3.3857590
## [11,] -6.2657601
## [12,] -8.7357613
## [13,] -8.0957581
## [14,] -9.3057610
## [15,] -9.4357582
## [16,] -8.5157601
## [17,] -7.8657585
## [18,] -7.3857590
## [19,] -6.5857597
## [20,] -8.4157616
## [21,] -8.2957588
## [22,] -9.7757584
## [23,] -11.2757584
## [24,] -14.2457596
## [25,] -14.8357597
## [26,] -14.8357597
## [27,] -13.5057579
## [28,] -14.5557610# Create DataFrames to store error metrics for both methods
metrics_columns <- c('MAPE', 'MAE', 'RMSE', 'MSE', 'R2')
direct_metrics <- data.frame(matrix(NA, nrow = length(horizons), ncol = length(metrics_columns)), row.names = horizons)
colnames(direct_metrics) <- metrics_columns
rolling_metrics <- data.frame(matrix(NA, nrow = length(horizons), ncol = length(metrics_columns)), row.names = horizons)
colnames(rolling_metrics) <- metrics_columns
# Fill DataFrames with the accuracy metrics from the results
for (horizon in horizons) {
direct_metrics[as.character(horizon), ] <- unlist(non_rolling_results$Results[[as.character(horizon)]]$Accuracy)
rolling_metrics[as.character(horizon), ] <- unlist(rolling_results$Results[[as.character(horizon)]]$Accuracy)
}
# Print the tables directly in R console
print(direct_metrics)## MAPE MAE RMSE MSE R2
## 7 0.08506050 4.942273 5.271315 27.78676 -7.1904761
## 14 0.04458097 2.545381 3.411227 11.63647 -0.3136622
## 21 0.14526581 7.959041 8.818567 77.76713 -4.1401832
## 28 0.12271716 6.674346 8.244068 67.96465 -1.8659661print(rolling_metrics)## MAPE MAE RMSE MSE R2
## 7 0.08516772 4.948978 5.280621 27.88496 -7.219420
## 14 0.04643121 2.624479 3.268669 10.68420 -0.206158
## 21 0.14738989 8.078537 8.966159 80.39201 -4.313680
## 28 0.12780238 6.934688 8.473741 71.80429 -2.027878knitr::kable(direct_metrics, caption = "Direct Forecast Error Metrics")| MAPE | MAE | RMSE | MSE | R2 | |
|---|---|---|---|---|---|
| 7 | 0.0850605 | 4.942273 | 5.271315 | 27.78676 | -7.1904761 |
| 14 | 0.0445810 | 2.545381 | 3.411227 | 11.63647 | -0.3136622 |
| 21 | 0.1452658 | 7.959041 | 8.818567 | 77.76713 | -4.1401832 |
| 28 | 0.1227172 | 6.674346 | 8.244068 | 67.96465 | -1.8659661 |
knitr::kable(rolling_metrics, caption = "Rolling Forecast Error Metrics")| MAPE | MAE | RMSE | MSE | R2 | |
|---|---|---|---|---|---|
| 7 | 0.0851677 | 4.948978 | 5.280621 | 27.88496 | -7.219420 |
| 14 | 0.0464312 | 2.624479 | 3.268669 | 10.68420 | -0.206158 |
| 21 | 0.1473899 | 8.078537 | 8.966159 | 80.39201 | -4.313680 |
| 28 | 0.1278024 | 6.934688 | 8.473741 | 71.80429 | -2.027878 |
close_prices <- df_TGLS$TGLS.Close
Se observa que los pronosticos no son muy acertados, a medida que
aumenta el horizonte la diferencia entre el pronóstico y el real
aumenta.
Basándonos en los resultados proporcionados, ninguno de los dos modelos (rolling o directo) parece ser significativamente mejor que el otro. Ambos muestran errores similares en términos de MAPE, MAE, RMSE y MSE para diferentes horizontes de predicción, y ambos tienen coeficientes de determinación (R2) negativos en todos los casos, lo que indica que ninguno de los modelos ajusta bien los datos.
Dado que ambos modelos muestran un desempeño deficiente en términos de ajuste a los datos, no podemos concluir que uno sea mejor que el otro en este momento.
##Best BIC
best_bic <- Inf # Mejor BIC inicializado a infinito
best_order_bic <- NULL # Mejor orden inicialmente a NULL
best_mdl_bic <- NULL # Mejor modelo inicialmente a NULL
pq_rng <- 0:7 # Rango para los parámetros p y q
d_rng <- 0:2 # Rango para el parámetro d
for (i in pq_rng) {
for (d in d_rng) {
for (j in pq_rng) {
tryCatch({
# Ajustar el modelo ARIMA con la orden (i, d, j)
tmp_mdl <- arima(train, order = c(i, d, j))
tmp_bic <- BIC(tmp_mdl) # Obtener el BIC del modelo temporal
# Verificar si el BIC actual es el mejor encontrado hasta ahora
if (tmp_bic < best_bic) {
best_bic <- tmp_bic # Actualizar el mejor BIC
best_order_bic <- c(i, d, j) # Guardar la mejor orden
best_mdl_bic <- tmp_mdl # Guardar el mejor modelo
}
}, error = function(e) {})
}
}
}## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1## Warning in log(s2): NaNs produced## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1## Warning in log(s2): NaNs produced
## Warning in log(s2): NaNs produced
## Warning in log(s2): NaNs produced## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1## Warning in log(s2): NaNs produced## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1best_order_bic## [1] 0 1 0best_order_bic <- c(0,1,0)TGLS_ts = ts(TGLS_xts, frequency = 1)
# Define los horizontes
horizons <- c(7, 14, 21, 28)
# Ejecutar el modelo ARIMA para cada horizonte y evaluar
non_rolling_results_bic <- non_rolling_arima_forecast(TGLS_ts, horizons, best_order_bic)
print(non_rolling_results_bic)## $Results
## $Results$`7`
## $Results$`7`$Forecast
## Time Series:
## Start = 2009
## End = 2015
## Frequency = 1
## [1] 52.46 52.46 52.46 52.46 52.46 52.46 52.46
##
## $Results$`7`$Accuracy
## $Results$`7`$Accuracy$MAPE
## [1] 0.08595714
##
## $Results$`7`$Accuracy$MAE
## [1] 4.994286
##
## $Results$`7`$Accuracy$RMSE
## [1] 5.323107
##
## $Results$`7`$Accuracy$MSE
## [1] 28.33546
##
## $Results$`7`$Accuracy$R2
## [1] -7.352212
##
##
##
## $Results$`14`
## $Results$`14`$Forecast
## Time Series:
## Start = 2002
## End = 2015
## Frequency = 1
## [1] 53.47 53.47 53.47 53.47 53.47 53.47 53.47 53.47 53.47 53.47 53.47 53.47
## [13] 53.47 53.47
##
## $Results$`14`$Accuracy
## $Results$`14`$Accuracy$MAPE
## [1] 0.04647383
##
## $Results$`14`$Accuracy$MAE
## [1] 2.627857
##
## $Results$`14`$Accuracy$RMSE
## [1] 3.278515
##
## $Results$`14`$Accuracy$MSE
## [1] 10.74866
##
## $Results$`14`$Accuracy$R2
## [1] -0.2134355
##
##
##
## $Results$`21`
## $Results$`21`$Forecast
## Time Series:
## Start = 1995
## End = 2015
## Frequency = 1
## [1] 45.46 45.46 45.46 45.46 45.46 45.46 45.46 45.46 45.46 45.46 45.46 45.46
## [13] 45.46 45.46 45.46 45.46 45.46 45.46 45.46 45.46 45.46
##
## $Results$`21`$Accuracy
## $Results$`21`$Accuracy$MAPE
## [1] 0.1396785
##
## $Results$`21`$Accuracy$MAE
## [1] 7.670477
##
## $Results$`21`$Accuracy$RMSE
## [1] 8.59692
##
## $Results$`21`$Accuracy$MSE
## [1] 73.90704
##
## $Results$`21`$Accuracy$R2
## [1] -3.885042
##
##
##
## $Results$`28`
## $Results$`28`$Forecast
## Time Series:
## Start = 1988
## End = 2015
## Frequency = 1
## [1] 44.25 44.25 44.25 44.25 44.25 44.25 44.25 44.25 44.25 44.25 44.25 44.25
## [13] 44.25 44.25 44.25 44.25 44.25 44.25 44.25 44.25 44.25 44.25 44.25 44.25
## [25] 44.25 44.25 44.25 44.25
##
## $Results$`28`$Accuracy
## $Results$`28`$Accuracy$MAPE
## [1] 0.1261087
##
## $Results$`28`$Accuracy$MAE
## [1] 6.848928
##
## $Results$`28`$Accuracy$RMSE
## [1] 8.403702
##
## $Results$`28`$Accuracy$MSE
## [1] 70.62221
##
## $Results$`28`$Accuracy$R2
## [1] -1.978031
##
##
##
##
## $Residuals
## $Residuals$`7`
## Time Series:
## Start = 2009
## End = 2015
## Frequency = 1
## [1] -1.480000 -2.980000 -5.950001 -6.540001 -6.540001 -5.209999 -6.260002
##
## $Residuals$`14`
## Time Series:
## Start = 2002
## End = 2015
## Frequency = 1
## [1] -0.1299973 0.7900009 1.4400024 1.9200020 2.7200012 0.8899994
## [7] 1.0100021 -0.4699974 -1.9699974 -4.9399986 -5.5299988 -5.5299988
## [13] -4.1999969 -5.2500000
##
## $Residuals$`21`
## Time Series:
## Start = 1995
## End = 2015
## Frequency = 1
## [1] -0.34000015 0.04000092 -2.09000015 -4.97000122 -7.44000244
## [6] -6.79999924 -8.01000214 -8.13999939 -7.22000122 -6.56999969
## [11] -6.09000015 -5.29000092 -7.12000275 -7.00000000 -8.47999954
## [16] -9.97999954 -12.95000076 -13.54000092 -13.54000092 -12.20999908
## [21] -13.26000214
##
## $Residuals$`28`
## Time Series:
## Start = 1988
## End = 2015
## Frequency = 1
## [1] -0.2799988 -0.7599983 -0.6300011 -0.4000015 -0.7500000 -1.3300018
## [7] -1.2099991 -1.5499992 -1.1699982 -3.2999992 -6.1800003 -8.6500015
## [13] -8.0099983 -9.2200012 -9.3499985 -8.4300003 -7.7799988 -7.2999992
## [19] -6.5000000 -8.3300018 -8.2099991 -9.6899986 -11.1899986 -14.1599998
## [25] -14.7500000 -14.7500000 -13.4199982 -14.4700012TGLS_ts = ts(TGLS_xts, frequency = 1)
# Define los horizontes
horizons <- c(7, 14, 21, 28)
# Ejecutar el modelo ARIMA para cada horizonte y evaluar
rolling_results_bic <- rolling_arima_forecast(TGLS_ts, horizons, best_order_bic)
rolling_results_bic## $Results
## $Results$`7`
## $Results$`7`$Forecast
## Time Series:
## Start = 2009
## End = 2015
## Frequency = 1
## [1] 52.46 52.46 52.46 52.46 52.46 52.46 52.46
##
## $Results$`7`$Accuracy
## $Results$`7`$Accuracy$MAPE
## [1] 0.08595714
##
## $Results$`7`$Accuracy$MAE
## [1] 4.994286
##
## $Results$`7`$Accuracy$RMSE
## [1] 5.323107
##
## $Results$`7`$Accuracy$MSE
## [1] 28.33546
##
## $Results$`7`$Accuracy$R2
## [1] -7.352212
##
##
##
## $Results$`14`
## $Results$`14`$Forecast
## Time Series:
## Start = 2002
## End = 2015
## Frequency = 1
## [1] 53.47 53.47 53.47 53.47 53.47 53.47 53.47 53.47 53.47 53.47 53.47 53.47
## [13] 53.47 53.47
##
## $Results$`14`$Accuracy
## $Results$`14`$Accuracy$MAPE
## [1] 0.04647383
##
## $Results$`14`$Accuracy$MAE
## [1] 2.627857
##
## $Results$`14`$Accuracy$RMSE
## [1] 3.278515
##
## $Results$`14`$Accuracy$MSE
## [1] 10.74866
##
## $Results$`14`$Accuracy$R2
## [1] -0.2134355
##
##
##
## $Results$`21`
## $Results$`21`$Forecast
## Time Series:
## Start = 1995
## End = 2015
## Frequency = 1
## [1] 45.46 45.46 45.46 45.46 45.46 45.46 45.46 45.46 45.46 45.46 45.46 45.46
## [13] 45.46 45.46 45.46 45.46 45.46 45.46 45.46 45.46 45.46
##
## $Results$`21`$Accuracy
## $Results$`21`$Accuracy$MAPE
## [1] 0.1396785
##
## $Results$`21`$Accuracy$MAE
## [1] 7.670477
##
## $Results$`21`$Accuracy$RMSE
## [1] 8.59692
##
## $Results$`21`$Accuracy$MSE
## [1] 73.90704
##
## $Results$`21`$Accuracy$R2
## [1] -3.885042
##
##
##
## $Results$`28`
## $Results$`28`$Forecast
## Time Series:
## Start = 1988
## End = 2015
## Frequency = 1
## [1] 44.25 44.25 44.25 44.25 44.25 44.25 44.25 44.25 44.25 44.25 44.25 44.25
## [13] 44.25 44.25 44.25 44.25 44.25 44.25 44.25 44.25 44.25 44.25 44.25 44.25
## [25] 44.25 44.25 44.25 44.25
##
## $Results$`28`$Accuracy
## $Results$`28`$Accuracy$MAPE
## [1] 0.1261087
##
## $Results$`28`$Accuracy$MAE
## [1] 6.848928
##
## $Results$`28`$Accuracy$RMSE
## [1] 8.403702
##
## $Results$`28`$Accuracy$MSE
## [1] 70.62221
##
## $Results$`28`$Accuracy$R2
## [1] -1.978031
##
##
##
##
## $Residuals
## $Residuals$`7`
## Time Series:
## Start = 2009
## End = 2015
## Frequency = 1
## test
## [1,] -1.480000
## [2,] -2.980000
## [3,] -5.950001
## [4,] -6.540001
## [5,] -6.540001
## [6,] -5.209999
## [7,] -6.260002
##
## $Residuals$`14`
## Time Series:
## Start = 2002
## End = 2015
## Frequency = 1
## test
## [1,] -0.1299973
## [2,] 0.7900009
## [3,] 1.4400024
## [4,] 1.9200020
## [5,] 2.7200012
## [6,] 0.8899994
## [7,] 1.0100021
## [8,] -0.4699974
## [9,] -1.9699974
## [10,] -4.9399986
## [11,] -5.5299988
## [12,] -5.5299988
## [13,] -4.1999969
## [14,] -5.2500000
##
## $Residuals$`21`
## Time Series:
## Start = 1995
## End = 2015
## Frequency = 1
## test
## [1,] -0.34000015
## [2,] 0.04000092
## [3,] -2.09000015
## [4,] -4.97000122
## [5,] -7.44000244
## [6,] -6.79999924
## [7,] -8.01000214
## [8,] -8.13999939
## [9,] -7.22000122
## [10,] -6.56999969
## [11,] -6.09000015
## [12,] -5.29000092
## [13,] -7.12000275
## [14,] -7.00000000
## [15,] -8.47999954
## [16,] -9.97999954
## [17,] -12.95000076
## [18,] -13.54000092
## [19,] -13.54000092
## [20,] -12.20999908
## [21,] -13.26000214
##
## $Residuals$`28`
## Time Series:
## Start = 1988
## End = 2015
## Frequency = 1
## test
## [1,] -0.2799988
## [2,] -0.7599983
## [3,] -0.6300011
## [4,] -0.4000015
## [5,] -0.7500000
## [6,] -1.3300018
## [7,] -1.2099991
## [8,] -1.5499992
## [9,] -1.1699982
## [10,] -3.2999992
## [11,] -6.1800003
## [12,] -8.6500015
## [13,] -8.0099983
## [14,] -9.2200012
## [15,] -9.3499985
## [16,] -8.4300003
## [17,] -7.7799988
## [18,] -7.2999992
## [19,] -6.5000000
## [20,] -8.3300018
## [21,] -8.2099991
## [22,] -9.6899986
## [23,] -11.1899986
## [24,] -14.1599998
## [25,] -14.7500000
## [26,] -14.7500000
## [27,] -13.4199982
## [28,] -14.4700012Best HQIC
best_hqic <- Inf # Mejor HQIC inicializado a infinito
best_order_hqic <- NULL # Mejor orden inicialmente a NULL
best_mdl_hqic <- NULL # Mejor modelo inicialmente a NULL
pq_rng <- 0:7 # Rango para los parámetros p y q
d_rng <- 0:2 # Rango para el parámetro d
for (i in pq_rng) {
for (d in d_rng) {
for (j in pq_rng) {
tryCatch({
# Ajustar el modelo ARIMA con la orden (i, d, j)
tmp_mdl <- arima(train, order = c(i, d, j))
tmp_hqic <- AIC(tmp_mdl, k = log(length(train))) # Obtener el HQIC del modelo temporal
# Verificar si el HQIC actual es el mejor encontrado hasta ahora
if (tmp_hqic < best_hqic) {
best_hqic <- tmp_hqic # Actualizar el mejor HQIC
best_order_hqic <- c(i, d, j) # Guardar la mejor orden
best_mdl_hqic <- tmp_mdl # Guardar el mejor modelo
}
}, error = function(e) {})
}
}
}## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1## Warning in log(s2): NaNs produced## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1## Warning in log(s2): NaNs produced
## Warning in log(s2): NaNs produced
## Warning in log(s2): NaNs produced## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1## Warning in log(s2): NaNs produced## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1
## Warning in arima(train, order = c(i, d, j)): possible convergence problem:
## optim gave code = 1TGLS_ts = ts(TGLS_xts, frequency = 1)
# Define los horizontes
horizons <- c(7, 14, 21, 28)
# Ejecutar el modelo ARIMA para cada horizonte y evaluar
non_rolling_results_hqic <- non_rolling_arima_forecast(TGLS_ts, horizons, best_order_hqic)TGLS_ts = ts(TGLS_xts, frequency = 1)
# Define los horizontes
horizons <- c(7, 14, 21, 28)
# Ejecutar el modelo ARIMA para cada horizonte y evaluar
rolling_results_hqic <- rolling_arima_forecast(TGLS_ts, horizons, best_order_hqic)Ejecutar pruebas para cada conjunto de residuales asociados a cada horizonte
fit_model <- arima(train, order=best_order)## Warning in arima(train, order = best_order): possible convergence problem:
## optim gave code = 1checkresiduals(fit_model)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(5,1,7)
## Q* = 11.342, df = 3, p-value = 0.01001
##
## Model df: 12. Total lags used: 15Los residuos se distribuyen alrededor de cero sin una tendencia clara. No se observa un patrón sistemático en los residuos a lo largo del tiempo. La ACF muestra que la mayoría de los residuos no están correlacionados con los rezagos anteriores. Solo hay algunas correlaciones significativas en los primeros rezagos. El histograma sugiere que los residuos siguen una distribución normal.
fit_model_bic <- arima(train, order=best_order_bic)
checkresiduals(fit_model_bic)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,0)
## Q* = 28.696, df = 10, p-value = 0.001395
##
## Model df: 0. Total lags used: 10Los residuos se distribuyen alrededor de cero sin una tendencia clara. No se observa un patrón sistemático en los residuos a lo largo del tiempo. La curva de distribución normal se ajusta bien a los datos. La ACF muestra que la mayoría de los residuos no están correlacionados con los rezagos anteriores.
fit_model_hqic <- arima(train, order=best_order_hqic)
checkresiduals(fit_model_hqic)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,0)
## Q* = 28.696, df = 10, p-value = 0.001395
##
## Model df: 0. Total lags used: 10El modelo ARIMA resultante de optimizar las métricas bic y hqic, es el mismo.