library(fpp3)
library(tictoc)

set.seed(334056)

dfLoadSar <- read.csv("C:\\Users\\Matteo\\OneDrive\\Desktop\\Tesi\\Dati\\LoadSar.csv")

date_parsed <- as.POSIXct(
  dfLoadSar$Date,
  format = "%Y-%m-%d %H:%M:%S",
  tz = "UTC"
)

date_errate <- dfLoadSar$Date[is.na(date_parsed)]
print(date_errate)

solo_data <- grepl("^\\d{4}-\\d{2}-\\d{2}$", dfLoadSar$Date)
dfLoadSar$Date[solo_data] <- paste0(dfLoadSar$Date[solo_data], " 00:00:00")

sum(is.na(dfLoadSar$Date))

dfLoadSar <- dfLoadSar[!is.na(dfLoadSar$Date), ]

sum(duplicated(dfLoadSar$Date))

rows_duplicates <- dfLoadSar %>%
  group_by(Date) %>%
  filter(n() > 1) %>%
  arrange(Date)

print(rows_duplicates)

dfLoadSar <- dfLoadSar %>%
  group_by(Date) %>%
  summarise(
    Total.Load..MW. = mean(Total.Load..MW., na.rm = TRUE),
    Forecast.Total.Load..MW. = mean(Forecast.Total.Load..MW., na.rm = TRUE),
    .groups = "drop"
  )

dfLoadSar$Date <- as.POSIXct(
  dfLoadSar$Date,
  format = "%Y-%m-%d %H:%M:%S",
  tz = "UTC"
)

LoadSarRaw <- as_tsibble(dfLoadSar, index = Date)

missingDate <- tsibble::scan_gaps(LoadSarRaw)
print(missingDate$Date)

LoadSarRaw <- LoadSarRaw |>
  fill_gaps() |>
  mutate(
    TotalLoad    = zoo::na.approx(Total.Load..MW.),
    ForecastLoad = zoo::na.approx(Forecast.Total.Load..MW.)
  )

LoadSarRaw <- LoadSarRaw %>%
  select(-Total.Load..MW., -Forecast.Total.Load..MW.)

lambdaG <- LoadSarRaw |>
  features(TotalLoad, guerrero) |>
  pull(lambda_guerrero)

LoadSar <- LoadSarRaw

initialStep <- 4 * 24 * 7 * 53
stepcv <- 4 * 24 * 8
teststepcv <- 4 * 24 * 8

cv <- LoadSar |>
  stretch_tsibble(
    .init = initialStep,
    .step = stepcv
  )

LoadSar |> autoplot(TotalLoad)

LoadSar |> gg_season(TotalLoad, period = "day") + ggtitle("day")
LoadSar |> gg_season(TotalLoad, period = "week") + ggtitle("Week")
LoadSar |> gg_season(TotalLoad, period = "year") + ggtitle("Year")

tic("fit")
cvEts <- cv |> model(ets = ETS(box_cox(TotalLoad, lambdaG)))
toc()

tic("forecast")
fcEtscv <- cvEts |>
  forecast(h = teststepcv)

fcEtscv |> accuracy(LoadSar)
toc()

tic("fit")
cvEtsStl <- cv |> model(
  etsStl = decomposition_model(
    STL(
      box_cox(TotalLoad, lambdaG) ~ season(period = "day") + season(period = "week") +
        season(period = "year"),
      robust = TRUE
    ),
    ETS(season_adjust ~ season("N"))
  )
)
toc()

tic("forecast")
fcEtsStlcv <- cvEtsStl |>
  forecast(h = teststepcv)

fcEtsStlcv |> accuracy(LoadSar)
toc()

tic("fit")
cvArima <- cv |> model(arima = ARIMA(box_cox(TotalLoad, lambdaG), stepwise = TRUE))
toc()

tic("forecast")
fcArimacv <- cvArima |>
  forecast(h = teststepcv)

fcArimacv |> accuracy(LoadSar)
toc()

tic("fit")
cvArimaStl <- cv |> model(
  arimaStl = decomposition_model(
    STL(
      box_cox(TotalLoad, lambdaG) ~ season(period = "day") + season(period = "week") +
        season(period = "year"),
      robust = TRUE
    ),
    ARIMA(season_adjust ~ PDQ(0, 0, 0) + pdq(3, 1, 3), stepwise = TRUE)
  )
)
toc()

tic("forecast")
fcArimaStlcv <- cvArimaStl |>
  forecast(h = teststepcv)

fcArimaStlcv |> accuracy(LoadSar)
toc()

horizonWeeks <- 24 * 8

etsRmse <- numeric(horizonWeeks)
etsStlRmse <- numeric(horizonWeeks)
arimaRmse <- numeric(horizonWeeks)
arimaStlRmse <- numeric(horizonWeeks)
prophetRmse <- numeric(horizonWeeks)
prophetStlRmse <- numeric(horizonWeeks)

for (i in 1:horizonWeeks) {
  teststepcv <- 4 * 24 * i
  fcEtscv <- cvEts |> forecast(h = teststepcv)
  etsRmse[i] <- fcEtscv |> accuracy(LoadSar) |> pull(RMSE)
}

for (i in 1:horizonWeeks) {
  teststep <- 4 * i
  fcEtsStlRidotto <- fcEtsStlcv |>
    group_by(.id, .model) |>
    filter(row_number() <= teststep) |>
    ungroup()
  
  etsStlRmse[i] <- fcEtsStlRidotto |> accuracy(LoadSar) |> pull(RMSE)
}

for (i in 1:horizonWeeks) {
  teststepcv <- 4 * 24 * i
  fcArimacv <- cvArima |> forecast(h = teststepcv)
  arimaRmse[i] <- fcArimacv |> accuracy(LoadSar) |> pull(RMSE)
}

for (i in 1:horizonWeeks) {
  teststepcv <- 4 * 24 * i
  fcArimaStlcv <- cvArimaStl |> forecast(h = teststepcv)
  arimaStlRmse[i] <- fcArimaStlcv |> accuracy(LoadSar) |> pull(RMSE)
}

for (i in 1:horizonWeeks) {
  teststep <- 4 * i
  fcArimaStlRidotto <- fcArimaStlcv |>
    group_by(.id, .model) |>
    filter(row_number() <= teststep) |>
    ungroup()
  
  arimaStlRmse[i] <- fcArimaStlRidotto |> accuracy(LoadSar) |> pull(RMSE)
}

horizonPlot <- 192

df_errori <- data.frame(
  Settimana = 1:horizonPlot,
  ARIMA_STL = arimaStlRmse[1:horizonPlot],
  ETS_STL = etsStlRmse[1:horizonPlot]
)

df_long <- df_errori |>
  pivot_longer(
    cols = c(ARIMA_STL, ETS_STL),
    names_to = "Modello",
    values_to = "RMSE"
  )

ggplot(data = df_long, aes(x = Settimana, y = RMSE, color = Modello)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  theme_minimal() +
  labs(
    title = "Confronto RMSE",
    subtitle = paste("Orizzonte di validazione:", horizonPlot, "Ore"),
    x = "Ore di forecast",
    y = "Errore (RMSE)"
  ) +
  scale_color_brewer(palette = "Set1") +
  scale_x_continuous(breaks = 1:horizonPlot)

dettaglio_cv_arima_stl <- fcArimaStlcv |>
  group_by(.id) |>
  accuracy(LoadSar) |>
  select(.id, .model, RMSE, MAE, MAPE)

print(dettaglio_cv_arima_stl, n = Inf)

dettaglio_cv_ets_stl <- fcEtsStlcv |>
  group_by(.id) |>
  accuracy(LoadSar) |>
  select(.id, .model, RMSE, MAE, MAPE)

print(dettaglio_cv_ets_stl, n = Inf)

library(stringr)

ids <- unique(cvArimaStl$.id)
modelli_estratti <- character(length(ids))

for (i in seq_along(ids)) {
  testo_report <- capture.output(
    cvArimaStl |> filter(.id == ids[i]) |> report()
  )
  
  riga_arima <- testo_report[str_detect(testo_report, "ARIMA\\(")]
  
  if (length(riga_arima) > 0) {
    modello_pulito <- str_replace(riga_arima[1], ".*Model:\\s*", "")
    modelli_estratti[i] <- str_trim(modello_pulito)
  } else {
    modelli_estratti[i] <- "Non trovato"
  }
}

tabella_frequenze <- data.frame(Modello_ARIMA = modelli_estratti) |>
  count(Modello_ARIMA, name = "Volte_Scelto") |>
  arrange(desc(Volte_Scelto))

print(tabella_frequenze)

ids <- unique(cvEtsStl$.id)
modelli_estratti_ets <- character(length(ids))

for (i in seq_along(ids)) {
  testo_report <- capture.output(
    cvEtsStl |> filter(.id == ids[i]) |> report()
  )
  
  riga_ets <- testo_report[str_detect(testo_report, "ETS\\(")]
  
  if (length(riga_ets) > 0) {
    modello_pulito <- str_replace(riga_ets[1], ".*Model:\\s*", "")
    modelli_estratti_ets[i] <- str_trim(modello_pulito)
  } else {
    modelli_estratti_ets[i] <- "Non trovato"
  }
}

tabella_frequenze_ets <- data.frame(Modello_ETS = modelli_estratti_ets) |>
  count(Modello_ETS, name = "Volte_Scelto") |>
  arrange(desc(Volte_Scelto))

print(tabella_frequenze_ets)

for (i in unique(cvArimaStl$.id)) {
  cat("\n==================================================\n")
  cat("  MODELLO SCELTO PER IL FOLD (id):", i, "\n")
  cat("==================================================\n")
  
  cvArimaStl |>
    filter(.id == i) |>
    report()
}