library(fpp3)

## US Leisure

leisure <- us_employment |>
  filter(Title == "Leisure and Hospitality", year(Month) > 2001) |>
  mutate(Employed = Employed / 1000) |>
  select(Month, Employed)

autoplot(leisure, Employed) +
  labs(y = "Employed (millions)")
fit <- leisure |>
  model(
    lm = TSLM(
      Employed ~
        trend(knot = yearmonth(c("2008 Jan", "2010 Jan"))) +
        fourier(K = 6)
    ),
    dhr = ARIMA(
      Employed ~
        trend(knot = yearmonth(c("2008 Jan", "2010 Jan"))) +
        fourier(K = 6)
    ),
    dhr2 = ARIMA(
      Employed ~
        trend(knot = yearmonth(c("2008 Jan", "2010 Jan"))) +
        fourier(K = 6) +
        pdq(0, 0, 0) +
        PDQ(0, 0, 0)
    )
  )
glance(fit)
fit |> select(dhr) |> report()

fit |>
  select(dhr) |>
  gg_tsresiduals()

fc <- fit |>
  forecast(h = 24)
autoplot(fc, leisure |> tail(12))


## 5-minute CALL CENTRE DATA ------------------------------------------------

calls <- readr::read_tsv("week12/callcenter.txt") |>
  rename(time = `...1`) |>
  pivot_longer(-time, names_to = "date", values_to = "volume") |>
  mutate(
    date = as.Date(date, format = "%d/%m/%Y"),
    datetime = as_datetime(date) + time
  ) |>
  as_tsibble(index = datetime)

calls |>
  fill_gaps() |>
  autoplot(volume)

calls |>
  fill_gaps() |>
  gg_season(volume, period = "day", alpha = 0.4) +
  guides(colour = "none")

calls_mdl <- calls |>
  mutate(idx = row_number()) |>
  update_tsibble(index = idx)
calls_fit <- calls_mdl |>
  model(ARIMA(volume ~ fourier(169, K = 10) + pdq(d = 0) + PDQ(0, 0, 0)))
report(calls_fit)

gg_tsresiduals(calls_fit, lag = 338)

calls_fit |>
  forecast(h = 1690) |>
  autoplot(calls_mdl |> filter(idx > 20000))

## US GASOLINE ---------------------------------------------------

us_gasoline |> autoplot(Barrels)

gasfit <- us_gasoline |>
  model(
    fourier1 = ARIMA(Barrels ~ fourier(K = 1) + PDQ(0, 0, 0)),
    fourier2 = ARIMA(Barrels ~ fourier(K = 2) + PDQ(0, 0, 0)),
    fourier3 = ARIMA(Barrels ~ fourier(K = 3) + PDQ(0, 0, 0)),
    fourier4 = ARIMA(Barrels ~ fourier(K = 4) + PDQ(0, 0, 0)),
    fourier5 = ARIMA(Barrels ~ fourier(K = 5) + PDQ(0, 0, 0)),
    fourier6 = ARIMA(Barrels ~ fourier(K = 6) + PDQ(0, 0, 0)),
    fourier7 = ARIMA(Barrels ~ fourier(K = 7) + PDQ(0, 0, 0)),
    fourier8 = ARIMA(Barrels ~ fourier(K = 8) + PDQ(0, 0, 0)),
    fourier9 = ARIMA(Barrels ~ fourier(K = 9) + PDQ(0, 0, 0)),
    fourier10 = ARIMA(Barrels ~ fourier(K = 10) + PDQ(0, 0, 0)),
    fourier11 = ARIMA(Barrels ~ fourier(K = 11) + PDQ(0, 0, 0)),
    fourier12 = ARIMA(Barrels ~ fourier(K = 12) + PDQ(0, 0, 0)),
    fourier13 = ARIMA(Barrels ~ fourier(K = 13) + PDQ(0, 0, 0)),
    fourier14 = ARIMA(Barrels ~ fourier(K = 14) + PDQ(0, 0, 0)),
    best_lm = TSLM(
      Barrels ~ trend(knots = yearweek(c("2006 W1", "2011 W1"))) +
        fourier(K = 6)
    ),
    best_lm2 = ARIMA(
      Barrels ~ trend(knots = yearweek(c("2006 W1", "2011 W1"))) +
        fourier(K = 6) +
        pdq(0, 0, 0) +
        PDQ(0, 0, 0)
    )
  )

glance(gasfit) |> arrange(AICc) |> select(.model, AICc)

gasfit |>
  select(fourier6) |>
  report()
gasfit |>
  select(fourier6) |>
  gg_tsresiduals()

gasfit |>
  select(fourier6) |>
  forecast(h = "3 years") |>
  autoplot(us_gasoline)

## Daily data with annual and weekly seasonality

vic_elec_daily <- vic_elec |>
  index_by(Date = date(Time)) |>
  summarise(
    Demand = sum(Demand) / 1e3,
    Temperature = max(Temperature),
    Holiday = any(Holiday)
  ) |>
  mutate(
    Day_Type = case_when(
      Holiday ~ "Holiday",
      wday(Date) %in% 2:6 ~ "Weekday",
      TRUE ~ "Weekend"
    )
  )

vic_elec_daily |>
  ggplot(aes(x = Temperature, y = Demand, colour = Day_Type)) +
  geom_point() +
  labs(x = "Maximum temperature", y = "Electricity demand (GW)")

vic_elec_daily |>
  pivot_longer(c(Demand, Temperature), names_to = "var", values_to = "value") |>
  ggplot(aes(x = Date, y = value)) +
  geom_line() +
  facet_grid(vars(var), scales = "free_y")

elec_fit <- vic_elec_daily |>
  model(
    ets = ETS(Demand),
    arima = ARIMA(log(Demand)),
    dhr = ARIMA(
      log(Demand) ~ Temperature +
        I(Temperature^2) +
        (Day_Type == "Weekday") +
        fourier(period = "year", K = 4)
    )
  )
elec_fit |>
  pivot_longer(ets:dhr, names_to = "model", values_to = "model_fit")

accuracy(elec_fit)

# ETS
elec_fit |>
  select(ets) |>
  report()
elec_fit |>
  select(ets) |>
  gg_tsresiduals()

elec_fit |>
  select(ets) |>
  augment() |>
  filter(Date <= "2014-03-31") |>
  ggplot(aes(x = Date, y = Demand)) +
  geom_line() +
  geom_line(aes(y = .fitted), col = "red")

# ARIMA
elec_fit |>
  select(arima) |>
  report()
elec_fit |>
  select(arima) |>
  gg_tsresiduals()

elec_fit |>
  select(arima) |>
  augment() |>
  filter(Date <= "2014-03-31") |>
  ggplot(aes(x = Date, y = Demand)) +
  geom_line() +
  geom_line(aes(y = .fitted), col = "red")

# DHR
elec_fit |>
  select(dhr) |>
  report()
elec_fit |>
  select(dhr) |>
  gg_tsresiduals()

elec_fit |>
  select(dhr) |>
  augment() |>
  filter(Date <= "2014-03-31") |>
  ggplot(aes(x = Date, y = Demand)) +
  geom_line() +
  geom_line(aes(y = .fitted), col = "red")

# Forecast 14 days ahead
vic_elec_future <- new_data(vic_elec_daily, 14) |>
  mutate(
    Temperature = c(rep(32, 7), rep(25, 7)),
    Holiday = c(TRUE, rep(FALSE, 13)),
    Day_Type = case_when(
      Holiday ~ "Holiday",
      wday(Date) %in% 2:6 ~ "Weekday",
      TRUE ~ "Weekend"
    )
  )
elec_fit |>
  forecast(new_data = vic_elec_future) |>
  autoplot(vic_elec_daily |> tail(14), level = 80) +
  labs(y = "Electricity demand (GW)")

# Forecast a year ahead using last year's temperatures
vic_elec_future <- new_data(vic_elec_daily, 365) |>
  mutate(
    Temperature = tail(vic_elec_daily$Temperature, 365),
    Holiday = Date %in%
      as.Date(c(
        "2015-01-01",
        "2015-01-26",
        "2015-03-09",
        "2015-04-03",
        "2015-04-06",
        "2015-04-25",
        "2015-06-08",
        "2015-10-02",
        "2015-11-03",
        "2015-12-25"
      )),
    Day_Type = case_when(
      Holiday ~ "Holiday",
      wday(Date) %in% 2:6 ~ "Weekday",
      TRUE ~ "Weekend"
    )
  )
forecast(elec_fit, new_data = vic_elec_future) |>
  filter(.model == "dhr") |>
  autoplot(vic_elec_daily |> tail(365), level = 80) +
  labs(y = "Electricity demand (GW)")