Harmonic Regression

library(tidyverse)
library(fpp3)
library(plotly)
library(patchwork)

vic_elec

Vamos a agregar los datos por hora:

vic_elec |> 
  mutate(fecha_hora = floor_date(Time, unit = "hour")) |> 
  relocate(fecha_hora)

wday(today())

[1] 3

vic_elec_hour <- vic_elec |> 
  index_by(time = ~ floor_date(., unit = "hour")) |> 
  summarise(
    demand = sum(Demand),
    mean_temp = mean(Temperature),
    max_temp = max(Temperature),
    holiday = any(Holiday)
  ) |> 
  mutate(
    day_type = case_when(
     holiday                ~ "holiday",
     wday(time) %in% c(2:6) ~ "weekday",
     TRUE                   ~ "weekend"
    )
  )
vic_elec_hour

p <- vic_elec_hour |> 
  autoplot(demand, color = "blue") 

ggplotly(p, dynamicTicks = TRUE) |> 
  rangeslider()

https://rpubs.com/pbenavides/harmonic_reg

vic_elec_hour |> 
  ggplot(aes(x = mean_temp, y = demand, 
             color = holiday)) +
  geom_point(alpha = 0.5)

vic_elec_hour |> 
  ggplot(aes(x = mean_temp, y = demand, 
             color = day_type)) +
  geom_point(alpha = 0.5)

vic_fit <- vic_elec_hour |> 
  model(
    harmonic1 = TSLM(demand ~ mean_temp + I(mean_temp^2) + day_type + fourier(period = "year", K = 4) + 
  fourier(period = "week", K = 3) +
  fourier(period = "day", K = 3)),
  harmonic2 = TSLM(demand ~ mean_temp + I(mean_temp^2) + day_type + fourier(period = "year", K = 15) + 
  fourier(period = "week", K = 10) +
  fourier(period = "day", K = 10)),
  harmonic3 = TSLM(demand ~ mean_temp + I(mean_temp^2) + day_type + fourier(period = "year", K = 6) + 
  fourier(period = "week", K = 5) +
  fourier(period = "day", K = 4))
  )

# report(vic_fit)

vic_aug <- vic_fit |> augment()

# vic_aug

accuracy(vic_fit) |> 
  select(.model, .type, MAPE, RMSE, MAE, MASE) |> 
  arrange(MAPE)

p <- vic_elec_hour |> 
  autoplot(demand) +
  geom_line(data = vic_aug, aes(x = time, y = .fitted),
            color = "firebrick") +
  facet_wrap(~ .model, ncol = 1)

ggplotly(p, dynamicTicks = TRUE, height = 700) |> 
  rangeslider()

vic_fit |> 
  select(harmonic2) |> 
  gg_tsresiduals()

vic_aug |> 
  filter(.model == "harmonic2") |> 
  features(.innov, unitroot_nsdiffs, period = 24)

vic_aug |> 
  filter(.model == "harmonic2") |> 
  features(.innov |> difference(24), unitroot_ndiffs)

p1 <- vic_aug |> 
  filter(.model == "harmonic2") |> 
  ACF(.innov |> difference(24)) |> 
  autoplot()
p2 <- vic_aug |> 
  filter(.model == "harmonic2") |> 
  PACF(.innov |> difference(24)) |> 
  autoplot()

p1 / p2

tictoc::tic()
# vic_fit <- vic_elec_hour |> 
#   model(
#   harmonic2 = ARIMA(demand ~ mean_temp + I(mean_temp^2) + day_type + fourier(period = "year", K = 6) + 
#   fourier(period = "week", K = 5) +
#   fourier(period = "day", K = 4) + pdq(0:2,0:1,0:2) + PDQ(0:1,0:1,0:2)
#   )
#   )
tictoc::toc()

0 sec elapsed