624_hw5
Haoming Chen
2026-03-08
library(fpp3)
library(tidyverse)8.1
#import data
pigs_vic <- aus_livestock %>%
filter(State == "Victoria", Animal == "Pigs")#a.
fit <- pigs_vic %>%
model(ETS(Count ~ error("A") + trend("N") + season("N")))
fc <- fit %>%
forecast(h = 4)
fc## # A fable: 4 x 6 [1M]
## # Key: Animal, State, .model [1]
## Animal State .model Month
## <fct> <fct> <chr> <mth>
## 1 Pigs Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + season(\"… 2019 Jan
## 2 Pigs Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + season(\"… 2019 Feb
## 3 Pigs Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + season(\"… 2019 Mar
## 4 Pigs Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + season(\"… 2019 Apr
## # ℹ 2 more variables: Count <dist>, .mean <dbl>#b.
resid_sd <- fit %>%
augment() %>%
summarise(s = sd(.resid, na.rm = TRUE)) %>%
pull(s)
first_fc <- fc %>%
as_tibble() %>%
slice(1) %>%
pull(.mean)
s <- fit %>%
glance() %>%
pull(sigma2) %>%
sqrt()
s## [1] 9353.115# my 95% interval
manual_interval <- tibble(
lower_95 = first_fc - 1.96 * s,
upper_95 = first_fc + 1.96 * s
)
manual_interval## # A tibble: 1 × 2
## lower_95 upper_95
## <dbl> <dbl>
## 1 76854. 113519.# R 95% interval
r_interval <- fc %>%
hilo(level = 95)
r_interval## # A tsibble: 4 x 7 [1M]
## # Key: Animal, State, .model [1]
## Animal State .model Month
## <fct> <fct> <chr> <mth>
## 1 Pigs Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + season(\"… 2019 Jan
## 2 Pigs Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + season(\"… 2019 Feb
## 3 Pigs Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + season(\"… 2019 Mar
## 4 Pigs Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + season(\"… 2019 Apr
## # ℹ 3 more variables: Count <dist>, .mean <dbl>, `95%` <hilo>8.5
# I choose Australia
country_data <- global_economy %>%
filter(Country == "Australia")#a.
country_data %>%
autoplot(Exports) +
labs(
title = "Annual Exports for Australia",
y = "Exports",
x = "Year"
)
#b. fit ETS and forcast
fit_ann <- country_data %>%
model(
ANN = ETS(Exports ~ error("A") + trend("N") + season("N"))
)
report(fit_ann)## Series: Exports
## Model: ETS(A,N,N)
## Smoothing parameters:
## alpha = 0.5659948
##
## Initial states:
## l[0]
## 12.98943
##
## sigma^2: 1.3621
##
## AIC AICc BIC
## 257.3943 257.8387 263.5756fc_ann <- fit_ann %>%
forecast(h = 10)
fc_ann %>%
autoplot(country_data) +
labs(
title = "ETS(A,N,N) Forecasts for Australia Exports",
y = "Exports",
x = "Year"
)
#c. RMSE
acc_ann <- fit_ann %>%
accuracy()
acc_ann## # A tibble: 1 × 11
## Country .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <fct> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Australia ANN Training 0.232 1.15 0.914 1.09 5.41 0.928 0.928 0.0125rmse_ann <- acc_ann %>%
pull(RMSE)
rmse_ann## [1] 1.146794#d. fit ETS and compare
fit_aan <- country_data %>%
model(
AAN = ETS(Exports ~ error("A") + trend("A") + season("N"))
)
report(fit_aan)## Series: Exports
## Model: ETS(A,A,N)
## Smoothing parameters:
## alpha = 0.4407571
## beta = 1e-04
##
## Initial states:
## l[0] b[0]
## 12.74993 0.1371607
##
## sigma^2: 1.3395
##
## AIC AICc BIC
## 258.3124 259.4662 268.6146fc_aan <- fit_aan %>%
forecast(h = 10)
fc_aan %>%
autoplot(country_data) +
labs(
title = "ETS(A,A,N) Forecasts for Australia Exports",
y = "Exports",
x = "Year"
)
acc_aan <- fit_aan %>%
accuracy()
acc_aan## # A tibble: 1 × 11
## Country .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <fct> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Australia AAN Training -7.46e-7 1.12 0.893 -0.387 5.39 0.907 0.904 0.109rmse_aan <- acc_aan %>%
pull(RMSE)
rmse_aan## [1] 1.116727#put training accuracy side by side
bind_rows(acc_ann, acc_aan)## # A tibble: 2 × 11
## Country .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <fct> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Australia ANN Training 2.32e-1 1.15 0.914 1.09 5.41 0.928 0.928 0.0125
## 2 Australia AAN Training -7.46e-7 1.12 0.893 -0.387 5.39 0.907 0.904 0.109#e compare forecast
bind_rows(
fc_ann %>% mutate(Model = "ETS(A,N,N)"),
fc_aan %>% mutate(Model = "ETS(A,A,N)")
) %>%
autoplot(country_data, level = NULL) +
labs(
title = "Comparison of Forecasts",
y = "Exports",
x = "Year"
)
#f.Manual 95% prediction
first_ann <- fc_ann %>%
as_tibble() %>%
slice(1) %>%
transmute(
Model = "ETS(A,N,N)",
Forecast = .mean,
Lower_95_manual = .mean - 1.96 * rmse_ann,
Upper_95_manual = .mean + 1.96 * rmse_ann
)first_aan <- fc_aan %>%
as_tibble() %>%
slice(1) %>%
transmute(
Model = "ETS(A,A,N)",
Forecast = .mean,
Lower_95_manual = .mean - 1.96 * rmse_aan,
Upper_95_manual = .mean + 1.96 * rmse_aan
)
manual_intervals <- bind_rows(first_ann, first_aan)
manual_intervals## # A tibble: 2 × 4
## Model Forecast Lower_95_manual Upper_95_manual
## <chr> <dbl> <dbl> <dbl>
## 1 ETS(A,N,N) 20.6 18.4 22.9
## 2 ETS(A,A,N) 20.8 18.6 23.0#R produced 95% interval
r_intervals_ann <- fc_ann %>%
hilo(level = 95) %>%
as_tibble() %>%
slice(1) %>%
mutate(Model = "ETS(A,N,N)")
r_intervals_aan <- fc_aan %>%
hilo(level = 95) %>%
as_tibble() %>%
slice(1) %>%
mutate(Model = "ETS(A,A,N)")
bind_rows(r_intervals_ann, r_intervals_aan)## # A tibble: 2 × 7
## Country .model Year
## <fct> <chr> <dbl>
## 1 Australia ANN 2018
## 2 Australia AAN 2018
## # ℹ 4 more variables: Exports <dist>, .mean <dbl>, `95%` <hilo>, Model <chr>8.6
# China GDP
china_gdp <- global_economy %>%
filter(Country == "China")
lambda_bc <- china_gdp %>%
features(GDP, guerrero) %>%
pull(lambda_guerrero)
lambda_bc## [1] -0.03446284fit_default <- china_gdp %>%
model(Default = ETS(GDP))
fit_damped <- china_gdp %>%
model(Damped = ETS(GDP ~ trend("Ad")))china_gdp_bc <- china_gdp %>%
mutate(GDP_bc = box_cox(GDP, lambda_bc))
fit_boxcox <- china_gdp_bc %>%
model(BoxCox = ETS(GDP_bc))
fit_boxcox_damped <- china_gdp_bc %>%
model(BoxCox_Damped = ETS(GDP_bc ~ trend("Ad")))# forecast
fc_default <- fit_default %>% forecast(h = 30)
fc_damped <- fit_damped %>% forecast(h = 30)
fc_boxcox <- fit_boxcox %>% forecast(h = 30)
fc_boxcox_damped <- fit_boxcox_damped %>% forecast(h = 30)
# plot and compare
autoplot(china_gdp, GDP) +
autolayer(fc_default, GDP, series = "Default", level = NULL) +
autolayer(fc_damped, GDP, series = "Damped", level = NULL) +
autolayer(fc_boxcox, GDP_bc, series = "BoxCox", level = NULL) +
autolayer(fc_boxcox_damped, GDP_bc, series = "BoxCox + Damped", level = NULL) +
labs(
title = "China GDP Forecasts with Different ETS Options",
y = "GDP",
x = "Year"
)
## 8.7
gas <- aus_production %>% select(Quarter, Gas)
# ETS model
fit1 <- gas %>% model(ETS(Gas))
fc1 <- fit1 %>% forecast(h = 32)
# Damped trend model
fit2 <- gas %>% model(ETS(Gas ~ trend("Ad")))
fc2 <- fit2 %>% forecast(h = 32)
# Plot forecasts
autoplot(gas, Gas) +
autolayer(fc1, series = "ETS") +
autolayer(fc2, series = "Damped ETS")
Multiplicative seasonality is needed because the seasonal
fluctuations increase as the level of gas production
increases.
Damped trend usually produces more realistic long-term forecasts
8.8
set.seed(12345678)
myseries <- aus_retail |>
filter(`Series ID` == sample(aus_retail$`Series ID`,1))#a
autoplot(myseries, Turnover)
#b
fit_hw <- myseries |>
model(HW = ETS(Turnover ~ error("M") + trend("A") + season("M")))
fit_hw_damped <- myseries |>
model(HW_damped = ETS(Turnover ~ error("M") + trend("Ad") + season("M")))#c
accuracy(fit_hw)## # A tibble: 1 × 12
## State Industry .model .type ME RMSE MAE MPE MAPE MASE RMSSE
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Northern T… Clothin… HW Trai… -0.0128 0.613 0.450 -0.469 5.15 0.513 0.529
## # ℹ 1 more variable: ACF1 <dbl>accuracy(fit_hw_damped)## # A tibble: 1 × 12
## State Industry .model .type ME RMSE MAE MPE MAPE MASE RMSSE
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Northern T… Clothin… HW_da… Trai… 0.0398 0.616 0.444 -0.0723 5.06 0.507 0.531
## # ℹ 1 more variable: ACF1 <dbl>#d
fit_hw |> gg_tsresiduals()
8.9
# train/test
train <- myseries |> filter(year(Month) < 2011)
test <- myseries |> filter(year(Month) >= 2011)lambda_bc <- train |>
features(Turnover, guerrero) |>
pull(lambda_guerrero)train_bc <- train |>
mutate(Turnover_bc = box_cox(Turnover, lambda_bc))
test_bc <- test |>
mutate(Turnover_bc = box_cox(Turnover, lambda_bc))fit_hw <- train |>
model(HW = ETS(Turnover ~ error("M") + trend("A") + season("M")))
fit_stl <- train_bc |>
model(
STL_ETS = decomposition_model(
STL(Turnover_bc),
ETS(season_adjust)
)
)fc_hw <- fit_hw |> forecast(new_data = test)
fc_stl <- fit_stl |> forecast(new_data = test_bc)
bind_rows(
accuracy(fc_hw, test),
accuracy(fc_stl, test_bc)
) |>
select(.model, RMSE)## # A tibble: 2 × 2
## .model RMSE
## <chr> <dbl>
## 1 HW 1.78
## 2 STL_ETS 0.405The STL + Box-Cox + ETS model performs much better because its RMSE is substantially lower. This means its forecasts are closer to the actual values on the test data.