Reentrega Time Series

  1. 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).

  2. 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.

  3. 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.

  4. 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(��,��).

  5. 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.

  6. 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.

Capturando la información del TICKER TGLS (Tecnoglass)

stock <- "TGLS" 
start_date <- as.Date("2014-04-01") 
end_date <- Sys.Date()

getSymbols(stock, src = "yahoo", 
           from = start_date, to = end_date)
[1] "TGLS"
TGLS_data <- na.omit(get(stock))  
TGLS<- data.frame(Date = index(TGLS_data), TGLS_data)  # Mostrar el marco de datos con las fechas usando datatable datatable(TGLS)
colnames(TGLS)
[1] "Date"          "TGLS.Open"     "TGLS.High"     "TGLS.Low"     
[5] "TGLS.Close"    "TGLS.Volume"   "TGLS.Adjusted"
plot(TGLS$Date, TGLS$TGLS.Close,       type = "l",       main = "Serie de Tiempo de TGLS.Close",      ylab = "Precio de Cierre",      xlab = "Fecha")  

hist(TGLS$TGLS.Close, main = "Histgram for Closing Price", xlab = "Freq", breaks = "Sturges", probability = TRUE) 
lines(density(TGLS$TGLS.Close))

La mayor parte del tiempo de la serie, el precio de cierre del stock ha estado alrededor de los 10 USD

TGLS$Month <- factor(month(TGLS$Date), labels = month.abb)  
ggplot(TGLS, aes(x = Month, y = TGLS.Close)) +   geom_boxplot() +   labs(x = "Mes", y = "Precio de Cierre", title = "Boxplot de TGLS.Close por Mes")  

De esta figura no se logra identificar una estacionariedad mensual, en general los datos parecen tener mayor dispersión por encima de la mediana. El mes de abril es en el que se encuentran la mayoría de datos atípicos.

Probando Estacionariedad con prueba de Dickey-Fuller

adf.test(Cl(TGLS_data$TGLS.Close), alternative = "stationary") 

    Augmented Dickey-Fuller Test

data:  Cl(TGLS_data$TGLS.Close)
Dickey-Fuller = -0.63585, Lag order = 13, p-value = 0.9757
alternative hypothesis: stationary

A partir de lo anterior, como P-value >0.05 podemos afirmar que esta serie de tiempo No es Estacionaria

Transformacion a serie estacionaria

* Usemos una diferenciación de primer orden:

price_diff <- diff(TGLS$TGLS.Close, lag = 1)

date_diff <- TGLS$Date[-1]  # Graficar las diferencias plot(date_diff, price_diff,       type = "l",       main = "Diferencia en Precio de Cierre de TGLS.Close",      ylab = "Diferencia de Precio de Cierre",      xlab = "Fecha") 

Realicemos una vez mas la prueba de Dickey Fuller para verificar si esta nueva serie es estacionaria

# Realiza la prueba ADF ignorando los valores faltantes 
adf.test(na.omit(price_diff), alternative = "stationary")  
Warning in adf.test(na.omit(price_diff), alternative = "stationary"): p-value
smaller than printed p-value

    Augmented Dickey-Fuller Test

data:  na.omit(price_diff)
Dickey-Fuller = -13.064, Lag order = 13, p-value = 0.01
alternative hypothesis: stationary

Como el P-valor < 0.05 entonces rechazamos la hipótesis nula, y aceptamos la hipótesis alternativa en la que la serie Si es estacionaria

price_diff<-na.omit(price_diff)

par(mfrow=c(1,2))
acf(price_diff, lag.max = 20) 
pacf(price_diff, lag.max = 20)

De estos graficos entendemos que una buena aproximación para los factores de diferenciación podrían ser AR(2) y diferenciación de nivel 1. (2,1,0)

TGLS_Price <- arima(price_diff, order = c(2, 1, 0))
summary(TGLS_Price)

Call:
arima(x = price_diff, order = c(2, 1, 0))

Coefficients:
          ar1      ar2
      -0.6397  -0.3780
s.e.   0.0184   0.0184

sigma^2 estimated as 0.5769:  log likelihood = -2898.81,  aic = 5803.63

Training set error measures:
                       ME      RMSE      MAE MPE MAPE      MASE        ACF1
Training set 0.0001253908 0.7593618 0.403262 NaN  Inf 0.8126558 -0.08701854
checkresiduals(TGLS_Price)


    Ljung-Box test

data:  Residuals from ARIMA(2,1,0)
Q* = 245.53, df = 8, p-value < 2.2e-16

Model df: 2.   Total lags used: 10

Este modelo no resulta muy bueno, dado que no pasa la prueba de Normalidad ni de Indepencia, no podemos afirmar que los residuales son ruido blanco, es decir, todavía existe un componente de autocorrelación que no se logró capturar en el modelo.

TGLS$Date <- as.Date(TGLS$Date)  # Crear objeto ts 
ts <- ts(TGLS$TGLS.Close, start = c(2014,04,01), end= c(2024,04,12), frequency = 365)
ts_info(ts)
 The ts series is a ts object with 1 variable and 3651 observations
 Frequency: 365 
 Start time: 2014 4 
 End time: 2024 4 
ts_decompose(ts)

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){})
            }
          }
        }
      }
    }
  }
best_aic 
[1] -6964.353
best_pdq 
[1] 2 1 2
best_PDQ
[1] 1 1 1

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)))
[1] "Train es:  3644"
test_w <- ts[train_size:(train_size + n_test-1)] 
print(paste("Test es: ", length(test_w)))
[1] "Test es:  7"
# 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
MAE MSE MAPE RMSE R2
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
MAE MSE MAPE RMSE R2
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)))
[1] "Train es:  3637"
test_w <- ts[train_size:(train_size + n_test-1)] 
print(paste("Test es: ", length(test_w)))
[1] "Test es:  14"
# 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
MAE MSE MAPE RMSE R2
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
MAE MSE MAPE RMSE R2
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)))
[1] "Train es:  3630"
test_w <- ts[train_size:(train_size + n_test-1)] 
print(paste("Test es: ", length(test_w)))
[1] "Test es:  21"
# 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
MAE MSE MAPE RMSE R2
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