library(readr)
library(fpp3)
library(zoo)
library(dplyr)
library(lubridate)
library(purrr)
library(rugarch)
library(tictoc)
library(stringr)

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

LoadSar <- LoadSarRaw |> rename(Load = TotalLoad)

LoadSar <- LoadSar %>%
  mutate(
    Giorno = date(Date),
    Ora_ = hour(Date),
    Quarto = as.integer(minute(Date) / 15),
    SlotQuarto = Ora_ * 4 + Quarto
  )

tsMakingLoadfpp3 <- function(dataset, ora) {
  
  dataset %>%
    as_tibble() %>%
    filter(SlotQuarto == ora) %>%
    arrange(Giorno) %>%
    mutate(
      diffLoad = Load - lag(Load)
    ) %>%
    filter(!is.na(diffLoad)) %>%
    as_tsibble(index = Giorno)
}

tsMakingLoad <- function(dataset, ora) {
  
  dataset %>%
    as_tibble() %>%
    filter(SlotQuarto == ora) %>%
    arrange(Giorno) %>%
    mutate(
      diffLoad = Load - lag(Load, 7)
    ) %>%
    filter(!is.na(diffLoad))
}

arimaLoadCV <- function(data_ts) {
  
  target_col <- "diffLoad"
  start_date <- "2023-01-01"
  refit_every <- 8
  
  start_date <- ymd(start_date)
  
  n_start <- data_ts %>%
    filter(Giorno < start_date) %>%
    nrow()
  
  cv_data <- data_ts %>%
    stretch_tsibble(
      .init = n_start,
      .step = refit_every
    )
  
  fit <- cv_data %>%
    model(
      arima = ARIMA(
        diffLoad ~ pdq(0, 0, 1) +
          PDQ(1, 0, 0, period = 7)
      )
    )
  
  fc <- fit %>%
    forecast(h = refit_every)
  
  fc_df <- fc %>%
    as_tibble() %>%
    select(
      .id,
      Giorno,
      Mu = .mean
    )
  
  actual_df <- data_ts %>%
    as_tibble() %>%
    select(
      Giorno,
      actual = all_of(target_col)
    )
  
  results <- fc_df %>%
    inner_join(actual_df, by = "Giorno") %>%
    mutate(
      errore = actual - Mu,
      abs_error = abs(errore),
      squared_error = errore^2
    ) %>%
    select(
      Giorno,
      actual,
      Mu,
      errore,
      abs_error,
      squared_error,
      .id
    )
  
  return(
    list(
      fit = fit,
      forecast = fc,
      results = results
    )
  )
}

etsLoadCV <- function(data_ts) {
  
  target_col <- "diffLoad"
  start_date <- "2023-01-01"
  refit_every <- 8
  
  start_date <- ymd(start_date)
  
  n_start <- data_ts %>%
    filter(Giorno < start_date) %>%
    nrow()
  
  cv_data <- data_ts %>%
    stretch_tsibble(
      .init = n_start,
      .step = refit_every
    )
  
  fit <- cv_data %>%
    model(
      ets = ETS(
        diffLoad ~ error("A") + trend("A") + season("A", period = 7)
      )
    )
  
  fc <- fit %>%
    forecast(h = refit_every)
  
  fc_df <- fc %>%
    as_tibble() %>%
    select(
      .id,
      Giorno,
      Mu = .mean
    )
  
  actual_df <- data_ts %>%
    as_tibble() %>%
    select(
      Giorno,
      actual = all_of(target_col)
    )
  
  results <- fc_df %>%
    inner_join(actual_df, by = "Giorno") %>%
    mutate(
      errore = actual - Mu,
      abs_error = abs(errore),
      squared_error = errore^2
    ) %>%
    select(
      Giorno,
      actual,
      Mu,
      errore,
      abs_error,
      squared_error,
      .id
    )
  
  return(
    list(
      fit = fit,
      forecast = fc,
      results = results
    )
  )
}

garchLoadCV <- function(data_ts) {
  refit_every <- 90
  target_col = "diffLoad"
  start_date = "2023-01-01"
  
  start_date <- ymd(start_date)
  
  n_start <- data_ts %>%
    filter(Giorno < start_date) %>%
    nrow()
  
  y <- as.numeric(data_ts[[target_col]])
  
  spec <- ugarchspec(
    variance.model = list(
      model = "sGARCH",
      garchOrder = c(1, 1)
    ),
    mean.model = list(
      armaOrder = c(0, 1),
      include.mean = TRUE
    ),
    distribution.model = "std"
  )
  
  roll <- ugarchroll(
    spec = spec,
    data = y,
    n.ahead = 1,
    n.start = n_start,
    refit.every = refit_every,
    refit.window = "recursive",
    solver = "hybrid",
    calculate.VaR = FALSE,
    keep.coef = TRUE
  )
  
  roll_df <- as.data.frame(roll)
  
  test_dates <- data_ts %>%
    slice((n_start + 1):n()) %>%
    select(Giorno, actual = all_of(target_col))
  
  results <- bind_cols(
    test_dates,
    roll_df
  )
  
  return(
    list(
      roll = roll,
      results = results
    )
  )
}

arima <- list()

tic("arima")
for (i in 0:95) {
  dataset <- tsMakingLoadfpp3(LoadSar, i)
  arima[[i + 1]] <- arimaLoadCV(dataset)
}
toc()

garch <- list()

tic("garch")
for (i in 0:95) {
  dataset <- tsMakingLoad(LoadSar, i)
  garch[[i + 1]] <- garchLoadCV(dataset)
}
toc()

ets <- list()

tic("ets")
for (i in 0:95) {
  dataset <- tsMakingLoadfpp3(LoadSar, i)
  ets[[i + 1]] <- etsLoadCV(dataset)
}
toc()

calcola_mae <- function(modello) {
  
  modello$results %>%
    as_tibble() %>%
    summarise(
      MAE = mean(abs(actual - Mu), na.rm = TRUE),
      MAPE = mean(abs((actual - Mu) / actual), na.rm = TRUE) * 100
    ) %>%
    pull(MAE)
}

crea_base_load <- function(dataset, slot, lag_base) {
  
  dataset %>%
    as_tibble() %>%
    filter(SlotQuarto == slot) %>%
    arrange(Giorno) %>%
    mutate(
      BaseLoad = lag(Load, lag_base)
    ) %>%
    select(
      Giorno,
      Load,
      BaseLoad
    )
}

calcola_mape_originale <- function(modello, dataset, slot, lag_base) {
  
  base_df <- crea_base_load(
    dataset = dataset,
    slot = slot,
    lag_base = lag_base
  )
  
  modello$results %>%
    as_tibble() %>%
    select(Giorno, Mu) %>%
    inner_join(base_df, by = "Giorno") %>%
    mutate(
      Load_hat = BaseLoad + Mu,
      errore_load = Load - Load_hat,
      APE = abs(errore_load / Load) * 100
    ) %>%
    summarise(
      MAPE = mean(APE, na.rm = TRUE)
    ) %>%
    pull(MAPE)
}

metriche_load_table <- tibble(
  SlotQuarto = 0:95,
  Ora_ = SlotQuarto %/% 4,
  Quarto = SlotQuarto %% 4,
  Orario = sprintf("%02d:%02d", Ora_, Quarto * 15),
  
  MAE_ARIMA = map_dbl(1:96, function(i) {
    calcola_mae(arima[[i]])
  }),
  
  MAE_GARCH = map_dbl(1:96, function(i) {
    calcola_mae(garch[[i]])
  }),
  
  MAE_ETS = map_dbl(1:96, function(i) {
    calcola_mae(ets[[i]])
  }),
  
  MAPE_ARIMA = map_dbl(1:96, function(i) {
    calcola_mape_originale(
      modello = arima[[i]],
      dataset = LoadSar,
      slot = i - 1,
      lag_base = 1
    )
  }),
  
  MAPE_GARCH = map_dbl(1:96, function(i) {
    calcola_mape_originale(
      modello = garch[[i]],
      dataset = LoadSar,
      slot = i - 1,
      lag_base = 7
    )
  }),
  
  MAPE_ETS = map_dbl(1:96, function(i) {
    calcola_mape_originale(
      modello = ets[[i]],
      dataset = LoadSar,
      slot = i - 1,
      lag_base = 1
    )
  })
)

print(metriche_load_table, n = Inf, width = Inf)

metriche_medie <- metriche_load_table %>%
  summarise(
    MAE_ARIMA_medio = mean(MAE_ARIMA, na.rm = TRUE),
    MAE_GARCH_medio = mean(MAE_GARCH, na.rm = TRUE),
    MAE_ETS_medio = mean(MAE_ETS, na.rm = TRUE),
    MAPE_ARIMA_medio = mean(MAPE_ARIMA, na.rm = TRUE),
    MAPE_GARCH_medio = mean(MAPE_GARCH, na.rm = TRUE),
    MAPE_ETS_medio = mean(MAPE_ETS, na.rm = TRUE)
  )

metriche_medie

mae_long <- metriche_load_table %>%
  select(SlotQuarto, Orario, MAE_ARIMA, MAE_GARCH, MAE_ETS) %>%
  pivot_longer(
    cols = c(MAE_ARIMA, MAE_GARCH, MAE_ETS),
    names_to = "Modello",
    values_to = "MAE"
  ) %>%
  mutate(
    Modello = recode(
      Modello,
      MAE_ARIMA = "ARIMA",
      MAE_GARCH = "GARCH",
      MAE_ETS = "ETS"
    )
  )

ggplot(mae_long, aes(x = SlotQuarto, y = MAE, color = Modello)) +
  geom_line(linewidth = 1) +
  scale_x_continuous(
    breaks = seq(0, 95, by = 4),
    labels = sprintf("%02d:00", 0:23)
  ) +
  labs(
    title = "MAE ARIMA vs GARCH vs ETS - Load Sardegna",
    x = "Time",
    y = "MAE",
    color = "Model"
  ) +
  theme_minimal()

mape_long <- metriche_load_table %>%
  select(SlotQuarto, Orario, MAPE_ARIMA, MAPE_GARCH, MAPE_ETS) %>%
  pivot_longer(
    cols = c(MAPE_ARIMA, MAPE_GARCH, MAPE_ETS),
    names_to = "Modello",
    values_to = "MAPE"
  ) %>%
  mutate(
    Modello = recode(
      Modello,
      MAPE_ARIMA = "ARIMA",
      MAPE_GARCH = "GARCH",
      MAPE_ETS = "ETS"
    )
  )

ggplot(mape_long, aes(x = SlotQuarto, y = MAPE, color = Modello)) +
  geom_line(linewidth = 1) +
  scale_x_continuous(
    breaks = seq(0, 95, by = 4),
    labels = sprintf("%02d:00", 0:23)
  ) +
  labs(
    title = "MAPE ARIMA vs GARCH vs ETS - Load Sardegna",
    x = "Time",
    y = "MAPE (%)",
    color = "Model"
  ) +
  theme_minimal()

slot_plot <- 48
data_inizio_plot <- as.Date("2023-01-01")
data_fine_plot <- as.Date("2023-12-31")

base_plot_ets <- crea_base_load(
  dataset = LoadSar,
  slot = slot_plot,
  lag_base = 1
)

plot_ets_df <- ets[[slot_plot + 1]]$results %>%
  as_tibble() %>%
  select(Giorno, Mu) %>%
  inner_join(base_plot_ets, by = "Giorno") %>%
  mutate(
    Load_ETS = BaseLoad + Mu
  ) %>%
  filter(
    Giorno >= data_inizio_plot,
    Giorno <= data_fine_plot
  ) %>%
  select(
    Giorno,
    Reale = Load,
    ETS = Load_ETS
  ) %>%
  pivot_longer(
    cols = c(Reale, ETS),
    names_to = "Serie",
    values_to = "Load"
  )

ggplot(plot_ets_df, aes(x = Giorno, y = Load, color = Serie)) +
  geom_line(linewidth = 1) +
  labs(
    title = paste0(
      "Real values vs ETS on LoadSar 2023 - Slot ",
      slot_plot,
      " (",
      sprintf("%02d:%02d", slot_plot %/% 4, (slot_plot %% 4) * 15),
      ")"
    ),
    x = "Date",
    y = "Load MW",
    color = "Graph"
  ) +
  theme_minimal()

crea_results_forecast_dataset <- function(dataset, slot, lag_base,
                                          data_inizio = NULL,
                                          data_fine = NULL) {
  
  results <- dataset %>%
    as_tibble() %>%
    filter(SlotQuarto == slot) %>%
    arrange(Giorno) %>%
    mutate(
      BaseLoad = lag(Load, lag_base),
      actual = Load - BaseLoad,
      Mu = ForecastLoad - BaseLoad,
      errore = actual - Mu,
      abs_error = abs(errore),
      squared_error = errore^2,
      errore_load = Load - ForecastLoad,
      abs_error_load = abs(errore_load),
      APE_load = abs(errore_load / Load) * 100
    ) %>%
    filter(
      !is.na(BaseLoad),
      !is.na(actual),
      !is.na(Mu),
      !is.na(Load),
      !is.na(ForecastLoad)
    )
  
  if (!is.null(data_inizio)) {
    results <- results %>%
      filter(Giorno >= as.Date(data_inizio))
  }
  
  if (!is.null(data_fine)) {
    results <- results %>%
      filter(Giorno <= as.Date(data_fine))
  }
  
  results %>%
    select(
      Giorno,
      Load,
      ForecastLoad,
      BaseLoad,
      actual,
      Mu,
      errore,
      abs_error,
      squared_error,
      errore_load,
      abs_error_load,
      APE_load
    )
}

calcola_metriche_forecast_dataset <- function(results_df) {
  
  results_df %>%
    summarise(
      MAE_diff = mean(abs_error, na.rm = TRUE),
      MAE_Load = mean(abs_error_load, na.rm = TRUE),
      MAPE_Load = mean(APE_load, na.rm = TRUE)
    )
}

calcola_metriche_terna <- function(dataset, slot, lag_base) {
  
  dataset %>%
    as_tibble() %>%
    filter(SlotQuarto == slot) %>%
    arrange(Giorno) %>%
    mutate(
      BaseLoad = lag(Load, lag_base),
      actual = Load - BaseLoad,
      Mu = ForecastLoad - BaseLoad,
      Load_hat = BaseLoad + Mu,
      errore_diff = actual - Mu,
      errore_load = Load - Load_hat,
      APE_load = abs(errore_load / Load) * 100
    ) %>%
    filter(
      !is.na(BaseLoad),
      !is.na(actual),
      !is.na(Mu),
      !is.na(Load),
      !is.na(ForecastLoad)
    ) %>%
    summarise(
      MAE_TERNA = mean(abs(errore_diff), na.rm = TRUE),
      MAPE_TERNA = mean(APE_load, na.rm = TRUE)
    )
}

metriche_load_table <- tibble(
  SlotQuarto = 0:95,
  Ora_ = SlotQuarto %/% 4,
  Quarto = SlotQuarto %% 4,
  Orario = sprintf("%02d:%02d", Ora_, Quarto * 15),
  
  MAE_ARIMA = map_dbl(1:96, function(i) {
    calcola_mae(arima[[i]])
  }),
  
  MAE_GARCH = map_dbl(1:96, function(i) {
    calcola_mae(garch[[i]])
  }),
  
  MAE_ETS = map_dbl(1:96, function(i) {
    calcola_mae(ets[[i]])
  }),
  
  MAE_TERNA = map_dbl(1:96, function(i) {
    calcola_metriche_terna(
      dataset = LoadSar,
      slot = i - 1,
      lag_base = 1
    ) %>%
      pull(MAE_TERNA)
  }),
  
  MAPE_ARIMA = map_dbl(1:96, function(i) {
    calcola_mape_originale(
      modello = arima[[i]],
      dataset = LoadSar,
      slot = i - 1,
      lag_base = 1
    )
  }),
  
  MAPE_GARCH = map_dbl(1:96, function(i) {
    calcola_mape_originale(
      modello = garch[[i]],
      dataset = LoadSar,
      slot = i - 1,
      lag_base = 7
    )
  }),
  
  MAPE_ETS = map_dbl(1:96, function(i) {
    calcola_mape_originale(
      modello = ets[[i]],
      dataset = LoadSar,
      slot = i - 1,
      lag_base = 1
    )
  }),
  
  MAPE_TERNA = map_dbl(1:96, function(i) {
    calcola_metriche_terna(
      dataset = LoadSar,
      slot = i - 1,
      lag_base = 1
    ) %>%
      pull(MAPE_TERNA)
  })
)

print(metriche_load_table, n = Inf, width = Inf)

tabella_medie_metriche <- metriche_load_table %>%
  summarise(
    MAE_TERNA = mean(MAE_TERNA, na.rm = TRUE),
    MAPE_TERNA = mean(MAPE_TERNA, na.rm = TRUE)
  )

print(tabella_medie_metriche, width = Inf)
