#load libraries
library(xgboost)
library(tidymodels)
## ── Attaching packages ────────────────────────────────────── tidymodels 1.2.0 ──
## ✔ broom 1.0.5 ✔ recipes 1.0.10
## ✔ dials 1.3.0 ✔ rsample 1.2.1
## ✔ dplyr 1.1.4 ✔ tibble 3.2.1
## ✔ ggplot2 3.5.1 ✔ tidyr 1.3.1
## ✔ infer 1.0.7 ✔ tune 1.2.1
## ✔ modeldata 1.4.0 ✔ workflows 1.1.4
## ✔ parsnip 1.2.1 ✔ workflowsets 1.1.0
## ✔ purrr 1.0.2 ✔ yardstick 1.3.1
## ── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
## ✖ purrr::discard() masks scales::discard()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::slice() masks xgboost::slice()
## ✖ recipes::step() masks stats::step()
## • Dig deeper into tidy modeling with R at https://www.tmwr.org
library(modeltime)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ readr 2.1.5
## ✔ lubridate 1.9.3 ✔ stringr 1.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ readr::col_factor() masks scales::col_factor()
## ✖ purrr::discard() masks scales::discard()
## ✖ dplyr::filter() masks stats::filter()
## ✖ stringr::fixed() masks recipes::fixed()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::slice() masks xgboost::slice()
## ✖ readr::spec() masks yardstick::spec()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(lubridate)
library(timetk)
library(fpp3)
## Registered S3 method overwritten by 'tsibble':
## method from
## as_tibble.grouped_df dplyr
## ── Attaching packages ──────────────────────────────────────────── fpp3 1.0.1 ──
## ✔ tsibble 1.1.5 ✔ feasts 0.4.1
## ✔ tsibbledata 0.4.1 ✔ fable 0.4.0
## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## ✖ fabletools::accuracy() masks yardstick::accuracy()
## ✖ lubridate::date() masks base::date()
## ✖ purrr::discard() masks scales::discard()
## ✖ dplyr::filter() masks stats::filter()
## ✖ fabletools::generate() masks infer::generate()
## ✖ fabletools::hypothesize() masks infer::hypothesize()
## ✖ tsibble::intersect() masks base::intersect()
## ✖ tsibble::interval() masks lubridate::interval()
## ✖ dplyr::lag() masks stats::lag()
## ✖ fabletools::null_model() masks parsnip::null_model()
## ✖ tsibble::setdiff() masks base::setdiff()
## ✖ dplyr::slice() masks xgboost::slice()
## ✖ tsibble::union() masks base::union()
library(patchwork)
data <- read.csv("multiTimeline new.csv", skip = 2, header = TRUE, fill = TRUE)
colnames(data) <- c('date','dodgers','ladodgers','otani')
data[data == "<1"] <- 0
new_data <- data |>
mutate(otani = as.numeric(data$otani),
date = yearmonth(date)) |>
as_tsibble(index = date)
la <- new_data |>
select(1,3) |>
mutate(
date = as.Date(date)
) |> as_tsibble(index = date)
otani <- new_data |>
select(1,4) |>
mutate(
date = as.Date(date)
)
otani |> autoplot()
## Plot variable not specified, automatically selected `.vars = otani`
#---------------------------------------------------------------
otani |> plot_time_series(
.date_var = date,
.value = otani
)
otani |> plot_seasonal_diagnostics(
.date_var = date,
.value = otani)
otani |> plot_anomaly_diagnostics(
.date_var = date,
.value = otani,
.alpha = 0.02
)
## frequency = 12 observations per 1 year
## trend = 60 observations per 5 years
otani |> plot_acf_diagnostics(
.date_var = date,
.value = otani)
## Max lag exceeds data available. Using max lag: 120
otani |> plot_stl_diagnostics(
.date_var = date,
.value = otani,
.facet_scales = 'free')
## frequency = 12 observations per 1 year
## trend = 60 observations per 5 years
#---------------------------------------------------------------
la |> plot_time_series(
.date_var = date,
.value = la$ladodgers
)
la |> plot_seasonal_diagnostics(
.date_var = date,
.value = la$ladodgers)
la |> plot_anomaly_diagnostics(
.date_var = date,
.value = la$ladodgers,
.alpha = 0.02
)
## frequency = 12 observations per 1 year
## trend = 60 observations per 5 years
la |> plot_acf_diagnostics(
.date_var = date,
.value = la$ladodgers)
## Max lag exceeds data available. Using max lag: 120
la |> plot_stl_diagnostics(
.date_var = date,
.value = la$ladodgers,
.facet_scales = 'free')
## frequency = 12 observations per 1 year
## trend = 60 observations per 5 years
#==============================================
long_data <- new_data %>%
pivot_longer(
cols = c(ladodgers, otani),
names_to = "source",
values_to = "value"
)
long_data |>
filter(date > as.Date("2020-10-31")) |>
ggplot(aes(x=date,y=value,color=source)) +
geom_point()+
geom_line()
## Warning: There was 1 warning in `filter()`.
## ℹ In argument: `date > as.Date("2020-10-31")`.
## Caused by warning:
## ! 不相容的方法 (">.vctrs_vctr", ">.Date") 用於 ">"
# *********************************************
#split the data into training and testing
la$date <- as.Date(la$date)
splits <- initial_time_split(la)
splits |>
tk_time_series_cv_plan() |>
plot_time_series_cv_plan(.date_var = date,
.value = ladodgers)
train <- training(splits)
test <- testing(splits)
## Create the (AUTO) models ----
### a. Auto ARIMA ----
arima_fit <- arima_reg() %>%
set_engine("auto_arima") %>%
fit(ladodgers ~ date, data = train)
## frequency = 12 observations per 1 year
### b. Boosted arima ----
arima_boost_fit <- arima_boost() %>%
set_engine("auto_arima_xgboost") %>%
fit(ladodgers ~ date, data = train)
## frequency = 12 observations per 1 year
### c. Exponential smoothing ----
ets_fit <- exp_smoothing() %>%
set_engine("ets") %>%
fit(ladodgers ~ date, data = train)
## frequency = 12 observations per 1 year
### d. Prophet ----
prophet_fit <- prophet_reg() %>%
set_engine("prophet") %>%
fit(ladodgers ~ date, data = train)
## Disabling weekly seasonality. Run prophet with weekly.seasonality=TRUE to override this.
## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
### e. Linear reg ----
lm_fit <- linear_reg() %>%
set_engine("lm") %>%
fit(ladodgers ~ date, data = train)
# *********************************************
## Add models to table ----
models_tbl <- modeltime_table(
arima_fit,
arima_boost_fit,
ets_fit,
prophet_fit,
lm_fit
)
# *********************************************
calibrate_tbl <- models_tbl %>%
modeltime_calibrate(new_data = test)
calibrate_tbl %>%
modeltime_forecast(
actual_data = la,
new_data = test
) %>%
plot_modeltime_forecast()
calibrate_tbl %>%
modeltime_accuracy()
## # A tibble: 5 × 9
## .model_id .model_desc .type mae mape mase smape rmse rsq
## <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 ARIMA(1,0,0)(0,1,1)[12] Test 8.65 50.0 1.58 42.4 15.3 0.323
## 2 2 ARIMA(1,0,0)(0,1,1)[12] Test 8.65 50.0 1.58 42.4 15.3 0.323
## 3 3 ETS(M,N,M) Test 10.0 63.8 1.84 56.5 15.7 0.293
## 4 4 PROPHET Test 8.11 78.7 1.48 48.9 11.7 0.564
## 5 5 LM Test 9.25 107. 1.69 54.7 17.2 0.235
refit_tbl <- calibrate_tbl %>%
modeltime_refit(data = la)
## frequency = 12 observations per 1 year
## frequency = 12 observations per 1 year
## frequency = 12 observations per 1 year
## Disabling weekly seasonality. Run prophet with weekly.seasonality=TRUE to override this.
## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
refit_tbl %>%
modeltime_forecast(h = "10 years", actual_data = la) %>%
plot_modeltime_forecast()
#========================================================
before <- la |>
filter(date < as.Date("2023-12-1")) |>
as_tsibble(index = date)
after <- la |>
filter(date >= as.Date("2023-12-1"))|>
as_tsibble(index = date)
before$date <- as.Date(before$date)
splits <- initial_time_split(before)
splits |>
tk_time_series_cv_plan() |>
plot_time_series_cv_plan(.date_var = date,
.value = ladodgers)
train <- training(splits)
test <- testing(splits)
## Create the (AUTO) models ----
### a. Auto ARIMA ----
arima_fit <- arima_reg() %>%
set_engine("auto_arima") %>%
fit(ladodgers ~ date, data = train)
## frequency = 12 observations per 1 year
### b. Boosted arima ----
arima_boost_fit <- arima_boost() %>%
set_engine("auto_arima_xgboost") %>%
fit(ladodgers ~ date, data = train)
## frequency = 12 observations per 1 year
### c. Exponential smoothing ----
ets_fit <- exp_smoothing() %>%
set_engine("ets") %>%
fit(ladodgers ~ date, data = train)
## frequency = 12 observations per 1 year
### d. Prophet ----
prophet_fit <- prophet_reg() %>%
set_engine("prophet") %>%
fit(ladodgers ~ date, data = train)
## Disabling weekly seasonality. Run prophet with weekly.seasonality=TRUE to override this.
## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
### e. Linear reg ----
lm_fit <- linear_reg() %>%
set_engine("lm") %>%
fit(ladodgers ~ date, data = train)
# *********************************************
## Add models to table ----
models_tbl <- modeltime_table(
arima_fit,
arima_boost_fit,
ets_fit,
prophet_fit,
lm_fit
)
# *********************************************
calibrate_tbl <- models_tbl %>%
modeltime_calibrate(new_data = test)
calibrate_tbl %>%
modeltime_forecast(
actual_data = la,
new_data = test
) %>%
plot_modeltime_forecast()
calibrate_tbl %>%
modeltime_accuracy()
## # A tibble: 5 × 9
## .model_id .model_desc .type mae mape mase smape rmse rsq
## <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 ARIMA(1,0,0)(0,1,1)[12] Test 4.99 48.7 0.806 33.1 10.2 0.447
## 2 2 ARIMA(1,0,0)(0,1,1)[12] Test 4.99 48.7 0.806 33.1 10.2 0.447
## 3 3 ETS(M,N,M) Test 5.91 81.6 0.956 41.3 13.0 0.404
## 4 4 PROPHET Test 6.03 99.3 0.974 52.7 8.57 0.588
## 5 5 LM Test 7.74 161. 1.25 64.2 12.2 0.0734
refit_tbl <- calibrate_tbl %>%
modeltime_refit(data = before)
## frequency = 12 observations per 1 year
## frequency = 12 observations per 1 year
## frequency = 12 observations per 1 year
## Disabling weekly seasonality. Run prophet with weekly.seasonality=TRUE to override this.
## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
refit_tbl %>%
modeltime_forecast(h = "1 years", actual_data = la) %>%
plot_modeltime_forecast()