library(fpp3)
library(fable.prophet)
library(tictoc)
library(ggplot2)
library(tidyr)
library(dplyr)

set.seed(334056)

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

dfPUN$DateTime <- as.POSIXct(
  paste(dfPUN$Data, sprintf("%02d:00:00", dfPUN$Ora - 1)),
  format = "%d/%m/%Y %H:%M:%S",
  tz = "UTC"
)

dfPUN <- dfPUN[, !names(dfPUN) %in% c("Data", "Ora")]

sum(is.na(dfPUN$DateTime))

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

date_errate <- dfPUN$DateTime[is.na(date_parsed)]
print(date_errate)

dfPUN$X..MWh <- as.numeric(gsub(",", ".", dfPUN$X..MWh))

sum(duplicated(dfPUN$DateTime))

dfPUN <- dfPUN %>%
  group_by(DateTime) %>%
  summarise(
    X..MWh = mean(X..MWh, na.rm = TRUE),
    .groups = "drop"
  )

dfPUN <- dfPUN |> rename(Price = X..MWh)

PUNRaw <- as_tsibble(dfPUN, index = DateTime)

tsibble::scan_gaps(PUNRaw)

PUNRaw <- PUNRaw |>
  fill_gaps() |>
  mutate(
    Price = zoo::na.approx(Price)
  )

lambdaG <- PUNRaw |>
  features(Price, guerrero) |>
  pull(lambda_guerrero)

PUN <- PUNRaw

initialStep <- 24 * 7 * 80
stepcv <- 24 * 7 * 80
teststepcv <- 24 * 7 * 1

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

PUNRaw |> autoplot(Price)

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

stl <- PUNRaw |> model(
  STL(
    Price ~ season(period = "day") + season(period = "week") +
      season(period = "year"),
    robust = TRUE
  )
)

stl |> components() |> autoplot()
stl |> gg_tsresiduals()

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

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

fcArimaStlcv |> accuracy(PUN)
toc()

cvArimaStl %>%
  as_tibble() %>%
  mutate(
    arima_specifico = map_chr(arimaStl, function(modello_singolo) {
      testo_report <- capture.output(report(modello_singolo))
      riga_arima <- testo_report[grepl("ARIMA\(", testo_report)]
      
      ```
      if (length(riga_arima) > 0) {
        return(str_extract(riga_arima[1], "ARIMA\\(.*"))
      } else {
        return("ARIMA non trovato")
      }
    })
    ```
    
  ) %>%
  count(arima_specifico, sort = TRUE)

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

print(dettaglio_cv_arima_stl)

trainRaw <- PUNRaw |> filter(DateTime < as.Date("2017-01-01"))
testRaw <- PUNRaw |> filter(DateTime >= as.Date("2017-01-01"))

train <- trainRaw
test <- testRaw

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

tic("forecast")
fcArimaStl <- fitArimaStl |>
  forecast(h = nrow(test))

fcArimaStl |> accuracy(PUN)
toc()

fcArimaStl |>
  autoplot(PUN)

fcArimaStl %>%
  filter(
    DateTime < as.POSIXct("2017-06-01")
  ) %>%
  autoplot(
    PUN %>%
      filter(
        DateTime >= as.POSIXct("2016-12-01"),
        DateTime < as.POSIXct("2017-06-01")
      )
  ) +
  labs(
    title = "ahead interval",
    y = "Price",
    x = "Date"
  )

horizonWeeks <- 7

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

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

plot(etsRmse, type = "l")

for (i in 1:horizonWeeks) {
  teststepcv <- 24 * i
  fcEtsStlcv <- cvEtsStl |> forecast(h = teststepcv)
  etsStlRmse[i] <- fcEtsStlcv |> accuracy(PUN) |> pull(RMSE)
}

plot(etsStlRmse, type = "l")

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

plot(arimaRmse, type = "l")

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

plot(arimaStlRmse, type = "l")

df_errori <- data.frame(
  Settimana = 1:horizonWeeks,
  ARIMA = arimaRmse,
  ARIMA_STL = arimaStlRmse,
  ETS = etsRmse,
  ETS_STL = etsStlRmse
)

df_long <- df_errori |>
  pivot_longer(
    cols = c(ARIMA, ARIMA_STL, ETS, 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: ARIMA vs ETS (con e senza STL)",
    subtitle = paste("Orizzonte di validazione:", horizonWeeks, "giorni"),
    x = "giorni di Forecast",
    y = "Errore (RMSE)"
  ) +
  scale_color_brewer(palette = "Set1") +
  scale_x_continuous(breaks = 1:horizonWeeks)

holidaysIt <- prophet::generated_holidays |>
  filter(country == "IT", year %in% 2004:2026) |>
  select(ds, holiday)

for (i in 1:horizonWeeks) {
  teststepcv <- 24 * 7 * i
  fcProphetcv <- cvProphet |> forecast(h = teststepcv)
  prophetRmse[i] <- fcProphetcv |> accuracy(PUN) |> pull(RMSE)
}

plot(prophetRmse, type = "l")

for (i in 1:horizonWeeks) {
  teststepcv <- 24 * 7 * i
  fcProphetStlcv <- cvProphetStl |> forecast(h = teststepcv)
  prophetStlRmse[i] <- fcProphetStlcv |> accuracy(PUN) |> pull(RMSE)
}

plot(prophetStlRmse, type = "l")

dettaglio_cv_arima <- fcArimacv |>
  group_by(.id) |>
  accuracy(PUN) |>
  select(.id, .model, RMSE, MAE, MAPE)

print(dettaglio_cv_arima)

dettaglio_cv_ets <- fcEtscv |>
  group_by(.id) |>
  accuracy(PUN) |>
  select(.id, .model, RMSE, MAE, MAPE)

print(dettaglio_cv_ets)

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

print(dettaglio_cv_ets_stl)

id_da_escludere <- c(12)

calcola_rmse_filtrato <- function(forecast_obj) {
  forecast_obj |>
    group_by(.id) |>
    accuracy(PUN) |>
    filter(!(.id %in% id_da_escludere)) |>
    summarise(RMSE_Medio = mean(RMSE, na.rm = TRUE))
}

rmse_ets <- calcola_rmse_filtrato(fcEtscv)
rmse_ets_stl <- calcola_rmse_filtrato(fcEtsStlcv)
rmse_arima <- calcola_rmse_filtrato(fcArimacv)
rmse_arima_stl <- calcola_rmse_filtrato(fcArimaStlcv)

tic("fit")
cvEtsStl <- cv |> model(
  etsStl = decomposition_model(
    STL(
      box_cox(Price, 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(PUN)
toc()

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

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

fcEtscv |> accuracy(PUN)
toc()

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

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

fcArimacv |> accuracy(PUN)
toc()

tic("fit")
cvProphet <- cv |> model(
  prophetModel = prophet(
    box_cox(Price, lambdaG) ~
      season(period = "day", order = 5, type = "multiplicative") +
      season(period = "week", order = 5, type = "multiplicative") +
      season(period = "year", order = 5, type = "multiplicative") +
      holiday(holidaysIt)
  )
)
toc()

tic("forecast")
fcProphetcv <- cvProphet |>
  forecast(h = 24 * 7 * 1)

fcProphetcv |> accuracy(PUN)
toc()

tic("fit")
cvProphetStl <- cv |> model(
  prophetStl = decomposition_model(
    STL(
      box_cox(Price, lambdaG) ~
        season(period = "day") +
        season(period = "week") +
        season(period = "year"),
      robust = TRUE
    ),
    prophet(season_adjust ~ holiday(holidaysIt))
  )
)
toc()

tic("forecast")
fcProphetStlcv <- cvProphetStl |>
  forecast(h = 24 * 7 * 1)

fcProphetStlcv |> accuracy(PUN)
toc()

tic("fit")
cvNnetar <- cv |> model(
  nnetar = NNETAR(Price ~ AR(p = 8, P = 1, period = "1 day"), n_nodes = 10)
)
toc()

tic("fit")
cvNnetarStl <- cv |> model(
  nnetarStl = decomposition_model(
    STL(
      Price ~ season(period = "day") + season(period = "week") +
        season(period = "year"),
      robust = TRUE
    ),
    NNETAR(season_adjust ~ AR(p = 8, P = 1, period = "1 day"), n_nodes = 10)
  )
)
toc()

tic("forecast")
fcNnetarcv <- cvNnetar |>
  forecast(h = teststepcv, simulate = FALSE)

fcNnetarcv |> accuracy(PUN)
toc()

tic("forecast")
fcNnetarStlcv <- cvNnetarStl |>
  forecast(h = teststepcv, simulate = FALSE)

fcNnetarStlcv |> accuracy(PUN)
toc()

tic("fit")
cvArimaStlHol <- cv |> model(
  arimaStlHol = decomposition_model(
    STL(
      TotalLoad ~ season(period = "day") + season(period = "week") +
        season(period = "year"),
      robust = TRUE
    ),
    ARIMA(season_adjust ~ PDQ(0, 0, 0) + pdq(d = 1) + isHoliday, stepwise = TRUE)
  )
)
toc()

tic("forecast")
fcArimaStlHolcv <- cvArimaStlHol |>
  forecast(new_data = testcv)

fcArimaStlHolcv |> accuracy(LoadItaly)
toc()

rmse

tic("fit")
cvArimaFourier <- cv |> model(
  arimaFourier = ARIMA(
    Price ~ PDQ(0, 0, 0) +
      fourier(period = "day", K = 5) +
      fourier(period = "week", K = 5) +
      fourier(period = "year", K = 5),
    stepwise = TRUE
  )
)
toc()

tic("forecast")
fcArimaFouriercv <- cvArimaFourier |>
  forecast(h = teststepcv)

fcArimaFouriercv |> accuracy(PUN)
toc()
