fable 包

library(tidyverse)
library(tsibble)
library(lubridate)
library(fable)
library(verification)
data(tourism)
#简化州名
tourism <- tourism %>%
  mutate(
    State = recode(State,
                   "Australian Capital Territory" = "ACT",
                   "New South Wales" = "NSW",
                   "Northern Territory" = "NT",
                   "Queensland" = "QLD",
                   "South Australia" = "SA",
                   "Tasmania" = "TAS",
                   "Victoria" = "VIC",
                   "Western Australia" = "WA"
    )
  )
#单序列预测
snowy <- tourism %>%
  filter(
    Region == "Snowy Mountains",
    Purpose == "Holiday"
  )
snowy %>%
  autoplot(Trips)

#seasonal naive、ETS和arima
fit <- snowy %>%
  model(
    snaive = SNAIVE(Trips ~ lag("year")),
    ets = ETS(Trips),
    arima = ARIMA(Trips)
  )

fit

fc <- fit %>%
  forecast(h = 12)

#该autoplot函数将产生所有预测的曲线图。默认情况下，level= C（80,95），因此80％和95％预测区间被示出。但是，为了避免混乱，我们将设置水平= NULL表明没有预测区间。
fc %>%
  autoplot(snowy, level = NULL) +
  ggtitle("Forecasts for Snowy Mountains holidays") +
  xlab("Year") +
  guides(colour = guide_legend(title = "Forecast"))

#如果要计算预测区间，可使用hilo函数
hilo(fc, level = 95)


#预测多系列
fit <- tourism %>%
  model(
    snaive = SNAIVE(Trips ~ lag("year")),
    ets = ETS(Trips),
    arima = ARIMA(Trips)
  )
fit

#选择特定的模型进行报告
fit %>%
  filter(Region == "Snowy Mountains", Purpose == "Holiday") %>%
  select(arima) %>%
  report()

fc <- fit %>%
  forecast(h = "3 years")
fc

fc %>%
  filter(Region == "Snowy Mountains") %>%
  autoplot(tourism, level = NULL) +
  xlab("Year") + ylab("Overnight trips (thousands)")

#模型预测准确率计算
train <- tourism %>%
  filter(year(Quarter) <= 2014)
#三个模型集成预测，简单计算平均值
fit <- train %>%
  model(
    ets = ETS(Trips),
    arima = ARIMA(Trips),
    snaive = SNAIVE(Trips)
  ) %>%
  mutate(mixed = (ets + arima + snaive) / 3)

fc <- fit %>% forecast(h = "3 years")

fc %>%
  filter(Region == "Snowy Mountains") %>%
  autoplot(tourism, level = NULL)

#预测与原始数据比较计算准确度
accuracy(fc, tourism)

#使用CRPS (Continuous Rank Probability Scores) and Winkler Scores (for 95% prediction intervals).计算准确性
fc_accuracy <- accuracy(fc, tourism,
                        measures = list(
                          point_accuracy_measures,
                          interval_accuracy_measures,
                          distribution_accuracy_measures
                        )
)

fc_accuracy %>%
  group_by(.model) %>%
  summarise(
    RMSE = mean(RMSE),
    MAE = mean(MAE),
    MASE = mean(MASE),
    Winkler = mean(winkler),
    CRPS = mean(CRPS)
  ) %>%
  arrange(RMSE)