## [1] "TGLS"

Comenzamos con analizar un poco los datos, previamente tratados para organizarlos correctamente segun una serie de tiempo. En el eje Y de esta grafica podemos ver la representacion de cada mes del año para su respectivo grafico de caja y bigote. Es facil de apreciar que todos los meses poseen una variabilidad parecida, con medianas y bigotes un poco diferentes.

 ts_plot((prices.ts), title =" Tecnoglass Adjusted Price",Ytitle= "price",Xtitle="Year")

con la funcion ts_plot graficamos la serie de tiempo

ts_plot(diff(prices.ts), title ="Diff Tecnoglass Adjusted Price",Ytitle= "price",Xtitle="Year")

En este lugar tenemos la grafica para la primera diferenciacion calculada sobre la serie. Vemos que fue util ya que logramos que los valores de la diferenciacion esten en torno al 0.

par(mfrow=c(1,2))
      acf(diff(prices.ts), lag.max = 20)
      pacf(diff(prices.ts), lag.max = 20)

Ahora tomamos la grafica de autocorrelacion y la grafica de autocorrelacion parcial. Vemos que en la primera grafica, a partir del segundo valor en x los valores en y se vuelven negativos. Luego, para la grafica de autocorrelacion parcial, se aprecian los 3 resagos, lo que se tomaria para el modelo ARIMA mas adelante.

best_ARIMA <- 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 = 12),
                             xreg=1:length(ts_in), 
                             method="CSS-ML")
                tmp_AIC <- AIC(fit) # Calculating HQIC instead of AIC
                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))
}

Con ayuda de esta funcion encontramos el mejor modelo ARIMA. Al correrlo para los distintos horizontes de tiempo se llega a lo siguiente: el mejor modelo arima es (3,1,1)(1,1,1).

 ####14
    splitted <- ts_split(prices.ts, sample.out = 14)
    train14<- splitted$train
    test14<-splitted$test
    mod14<- arima(train14, order = c(3,1,1), seasonal = list(order=c(1,1,1), period =12))
    mod14
## 
## Call:
## arima(x = train14, order = c(3, 1, 1), seasonal = list(order = c(1, 1, 1), period = 12))
## 
## Coefficients:
##           ar1      ar2      ar3     ma1     sar1     sma1
##       -0.6661  -0.1559  -0.4176  0.5139  -0.3875  -0.6080
## s.e.   0.1634   0.1107   0.0988  0.1689   0.1569   0.1346
## 
## sigma^2 estimated as 0.0005332:  log likelihood = 266.22,  aic = -518.43
    fcst14 <- forecast(mod14, h=14)
    ####21
    splitted <- ts_split(prices.ts, sample.out = 21)
    train21<- splitted$train
    test21<-splitted$test
    mod21<- arima(train21, order = c(3,1,1), seasonal = list(order=c(1,1,1), period =12))
    mod21
## 
## Call:
## arima(x = train21, order = c(3, 1, 1), seasonal = list(order = c(1, 1, 1), period = 12))
## 
## Coefficients:
##           ar1      ar2      ar3     ma1     sar1     sma1
##       -0.6772  -0.1696  -0.4417  0.5259  -0.3717  -0.5994
## s.e.   0.1510   0.1171   0.1033  0.1560   0.1628   0.1416
## 
## sigma^2 estimated as 0.0005385:  log likelihood = 249.5,  aic = -484.99
    fcst21 <- forecast(mod21, h=21)
    ####28
    splitted <- ts_split(prices.ts, sample.out = 28)
    train28<- splitted$train
    test28<-splitted$test
    mod28<- arima(train28, order = c(3,1,1), seasonal = list(order=c(1,1,1), period =12))
    mod28
## 
## Call:
## arima(x = train28, order = c(3, 1, 1), seasonal = list(order = c(1, 1, 1), period = 12))
## 
## Coefficients:
##           ar1      ar2      ar3     ma1     sar1     sma1
##       -0.7003  -0.1629  -0.4134  0.5474  -0.3310  -0.5969
## s.e.   0.1545   0.1190   0.1049  0.1607   0.1801   0.1660
## 
## sigma^2 estimated as 0.0005336:  log likelihood = 234.04,  aic = -454.07
    fcst28 <- forecast(mod28, h=28)

Luego, en esta parte hacemos la division de los datos en: train y test, correspondientes a la longitude del horizonte de cada caso, con ayuda de ‘ts_split’. Ademas, Generamos el modelo, teniendo en cuenta el mejor AIC, con las metricas halladas con ayuda de la funcion anterior: (3,1,1)(1,1,1).

p<- test_forecast((prices.ts), forecast.obj = fcst14, test = test14)
      p <- ggplotly(p, tooltip = "text", dynamicTicks = TRUE) %>% layout(title = "Forecast test 14 days")
      p
p <- test_forecast((prices.ts), forecast.obj = fcst21, test = test21)
      p <- ggplotly(p, tooltip = "text", dynamicTicks = TRUE) %>% layout(title = "Forecast test 21 days")
      p
p<- test_forecast((prices.ts), forecast.obj = fcst28, test = test28)
      p <- ggplotly(p, tooltip = "text", dynamicTicks = TRUE) %>% layout(title = "Forecast test 28 days")
      p

Al comparar las predicciones con los datos reales, es evidente que tienden a acercarse. Aunque muestran diferencias en su comportamiento, generalmente siguen una dirección similar. Esta convergencia entre las predicciones y los datos reales sugiere que el modelo utilizado para hacer las predicciones captura de alguna manera la dinámica subyacente del fenómeno en estudio

plot_forecast(fcst14,
                    title =" Tecnoglass stocks -  Forecast 14 days",
                    Ytitle = "Adjusted.Price",
                    Xtitle = "Year")
 plot_forecast(fcst21,
                    title =" Tecnoglass stocks -  Forecast 21 days",
                    Ytitle = " Adjusted.Price",
                    Xtitle = "Year")
plot_forecast(fcst28,
                    title =" Tecnoglass stocks -  Forecast 28 days",
                    Ytitle = " Adjusted.Price",
                    Xtitle = "Year")

Estas son las predicciones para los diferentes horizontes temporales del precio ajustado de Tecnoglass. Encontramos esos valores esperados segun los modelos obtenidos anteriormente.

checkresiduals(mod14)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(3,1,1)(1,1,1)[12]
## Q* = 22.739, df = 18, p-value = 0.2008
## 
## Model df: 6.   Total lags used: 24

Notemos que segun el orden de diferenciacion se logra que la serie gire al rededor del 0. Para la grafica de autocorrelacion se aprecia como no hay rezagos que lleguen a la banda de significancia.Finalmente, Los datos parecen ajustarse a una distribucion normal.

checkresiduals(mod21)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(3,1,1)(1,1,1)[12]
## Q* = 21.809, df = 18, p-value = 0.2406
## 
## Model df: 6.   Total lags used: 24

Para un horizonte de 21 dias sucede lo mismo con el orden de diferenciacion, logramos que a lo largo del tiempo los valores giren en torno al 0. Luego, para la autocorrelacion la grafica sigue sin tener rezagos. Adicionalmente, se mantiene el ajuste a una distibucion normal de los residuos.

checkresiduals(mod28)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(3,1,1)(1,1,1)[12]
## Q* = 21.267, df = 17, p-value = 0.2146
## 
## Model df: 6.   Total lags used: 23

Para el horizonte de los 28 dias tenemos residuos casi iguales a los anteriores, en cuanto a la distribucion de los datos y el comportamiento de las serie segun el grado de diferenciacion. No obstante, una pequeña diferencia es que se logra apreciar un rezago en el grafico de correlacion.

Comprobamos los residuos de los modelos. Obtenemos resultados favorables, ya que gracias al test de Ljung-Box se tiene lo siguienteç: para el modelo con un horizonte de 14 dias un p-value(0.208)>0.05, para el modelo con horizonte de 21 dias un p-value(0.240)>0.05 y, finalmente, para el modelo con horizonte de 28 dias un p-value(0.214)>0.05. Por tanto, comprobamos que para todos los modelos los residuos son ruido blanco, lo que nos corrobora los resultados.

mod_no<-auto.arima(prices.ts)
no_rolling_14<-forecast(mod_no, h=14)
      plot(no_rolling_14, main= "Forecast with horizon = 14")

 no_rolling_14<-forecast(mod_no, h=21)
      plot(no_rolling_14, main= "Forecast with horizon = 21")

 no_rolling_14<-forecast(mod_no, h=28)
      plot(no_rolling_14, main= "Forecast with horizon = 28")

Como podemos apreciar, para las tres predicciones, con sus respectivos horizontes, los valores futuros llegan a ser muy constantes, cosa que contrasta con la prediccion usando rolling.

res_no_rolling<- checkresiduals(mod_no)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(3,1,2)(1,0,0)[12]
## Q* = 18.808, df = 18, p-value = 0.4037
## 
## Model df: 6.   Total lags used: 24
print(res_no_rolling)
## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(3,1,2)(1,0,0)[12]
## Q* = 18.808, df = 18, p-value = 0.4037

Ahora vemos que la seria logra moverse con una media aproximada a 0. Se encuentra un rezago toca la banda de significancia. Y se aprecia la distribucion de los datos del modelo, un poco parecido a la que resulta del rolling. Luego, al realizar el Ljung Box-test se encuentra un p-valor(0.404) mayor al nivel de signficancia, por tanto, se corrobora que no hay evidencia significativa para afirmar que existe correlacion del modelo con los datos reales.

cat("Métricas para pronóstico con rolling a 14 días:\n")
## Métricas para pronóstico con rolling a 14 días:
      print(accuracy(fcst14))
##                       ME       RMSE        MAE        MPE      MAPE      MASE
## Training set 0.004568062 0.02203566 0.01378589 0.06004662 0.1811544 0.4134427
##                    ACF1
## Training set -0.0436084
      cat("\n")
      cat("Métricas para pronóstico con rolling a 21 días:\n")
## Métricas para pronóstico con rolling a 21 días:
      print(accuracy(fcst21))
##                       ME       RMSE        MAE        MPE      MAPE      MASE
## Training set 0.005053515 0.02207984 0.01369129 0.06639231 0.1801145 0.4351714
##                     ACF1
## Training set -0.05261843
      cat("\n")
      cat("Métricas para pronóstico con rolling a 28 días:\n")
## Métricas para pronóstico con rolling a 28 días:
      print(accuracy(fcst28))
##                      ME       RMSE        MAE       MPE      MAPE      MASE
## Training set 0.00491326 0.02191132 0.01343822 0.0646604 0.1770409 0.4199909
##                     ACF1
## Training set -0.06039411
cat("Métricas para pronósticos realizados sin rolling :\n")
## Métricas para pronósticos realizados sin rolling :
print(accuracy(no_rolling_14))
##                        ME       RMSE        MAE         MPE      MAPE      MASE
## Training set 0.0005531208 0.02051171 0.01048871 0.006659285 0.1373236 0.3187793
##                    ACF1
## Training set 0.01868474

Las métricas de rendimiento proporcionan una visión detallada del comportamiento del modelo de pronóstico en diferentes escenarios temporales. A pesar de las ligeras variaciones observadas en algunas métricas, como el MAPE, el modelo exhibe una estabilidad general en su capacidad para predecir los valores futuros. La consistencia en el MASE sugiere que el modelo supera consistentemente las expectativas de un enfoque de pronóstico simple. No obstante, la autocorrelación negativa del error sugiere que el modelo podría estar subestimando levemente la estructura de dependencia temporal en los datos. En resumen, las metricas indican que el modelo presenta un desempeño sólido y confiable en la tarea de pronóstico, pero se podrían explorar estrategias para mejorar su capacidad para capturar la correlación serial de los errores.

 par(mfrow=c(2,2))
      last_14_real <- tail(prices.ts, 14)
      plot(fcst14$mean, main ="forecast with horizon = 14 vs real data")
      points(last_14_real, col = "red")  
      
      last_21_real <- tail(prices.ts, 21)
      plot(fcst21$mean, main ="forecast with horizon = 21 vs real data")
      points(last_21_real, col = "red") 
      
      last_28_real <- tail(prices.ts, 28)
      plot(fcst28$mean, main ="forecast with horizon = 28 vs real data")
      points(last_28_real, col = "red")

Ahora vemos mas de cerca esa grafica de dispersion los valores reales reales comparados con las predicciones, en este caso usando rolling para las predicciones. las predicciones parecen seguir el mismo camino de los datos reales, aunque tome valores algo apartes.

Para finalizar, se intentó generar distintos modelos teniendo en cuenta otros parámetros como el HQIC; sin embargo, R no pudo encontrarlos.