## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
options(tibble.width = Inf, knitr.kable.NA = "")

# Helper to render fabletools tables as HTML-safe kable
kable <- function(x, ...) {
  x |>
    dplyr::mutate(dplyr::across(
      where(~ inherits(.x, "agg_vec")),
      ~ gsub("<", "&lt;", gsub(">", "&gt;", as.character(.x)))
    )) |>
    knitr::kable(...)
}

## ----libraries----------------------------------------------------------------
library(bayesRecon)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)
library(fable.bayesRecon)
library(ggplot2)
library(scales)

## ----data---------------------------------------------------------------------
data <- swiss_tourism$ts[, -1, drop = FALSE] |>             # Remove 1st column
  as_tsibble() |>                                           # Convert in tsibble
  rename(Month = index, Canton = key, Tourists = value) |>  # Rename
  aggregate_key(Canton, Tourists = sum(Tourists))           # Aggregate all cantons

data |>  head(5) |> kable()

## ----fit----------------------------------------------------------------------
fit <- data |>
  filter(Month < yearmonth("2024 Feb")) |>                  # Remove test set
  model(base = ETS(Tourists ~ trend("A") + season("A"))) |> # fit ETS model
  reconcile(t = bayesRecon_t(base),                         # Reconcile with t-Rec
            mint = min_trace(base))                         # Reconcile with MinT

fit |> head(5) |> kable()

## ----forecast-----------------------------------------------------------------
fc <- fit |>
  forecast(h = "1 year")

# Showing only the one step ahead forecast
fc |> filter(Month == yearmonth("Feb 2024")) |> head(5) |> kable()

## ----accuracy-----------------------------------------------------------------
results <- fc |>
  accuracy(data, measures = list(RMSSE = RMSSE, MASE = MASE, CRPS = CRPS)) |>
  group_by(.model) |>
  summarise(RMSSE = mean(RMSSE), MASE  = mean(MASE), CRPS  = mean(CRPS))

results |> kable(digits = 3)

## ----plot-aggregated, fig.width=7, fig.height=4-------------------------------
# Aggregated series
fc |>
  filter(is_aggregated(Canton)) |>
  autoplot(data, level=95) +
  scale_y_continuous(labels = label_number(suffix = "K", scale = 1e-3))+
  labs(title = "Aggregated Swiss Tourism", y = "Tourists", x = NULL) +
  theme_minimal()

## ----plot-aggregated-zoom, fig.width=7, fig.height=4--------------------------
# Zoom on the forecasts for the aggregated series
fc |>
  filter(is_aggregated(Canton)) |>
  autoplot(data |> filter(Month >= yearmonth("2022 Jan")), level = 95) +
  scale_y_continuous(labels = label_number(suffix = "K", scale = 1e-3))+
  labs(title = "Aggregated Swiss Tourism (Jan 2022 - Jan 2025)", y = "Tourists", x = NULL) +
  theme_minimal()

## ----plot-canton1, fig.width=7, fig.height=4----------------------------------
# Canton Ticino
fc |>
  filter(Canton == "Ticino") |>
  autoplot(data |> filter(Month >= yearmonth("2022 Jan")),
           level=95) +
  scale_y_continuous(labels = label_number(suffix = "K", scale = 1e-3))+
  labs(title = paste("Canton: Ticino"), y = "Tourists", x = NULL) +
  theme_minimal()

## ----plot-canton2, fig.width=7, fig.height=4----------------------------------
# Canton Luzern
fc |>
  filter(Canton == "Luzern") |>
  autoplot(data |> filter(Month >= yearmonth("2022 Jan"),
                                 Canton == "Luzern"),
           level = 80) +
  scale_y_continuous(labels = label_number(suffix = "K", scale = 1e-3))+
  labs(title = paste("Canton: Luzern"), y = "Tourists", x = NULL) +
  theme_minimal()