Forecast #4
Hunter Patenaude
Nov 4
Load packages and data
library(fpp3)
library(readxl)
# install and load any package necessaryData, Plot, and Decomp
Forecast_Bacon_Price <- read_excel("~/Downloads/Forecast Bacon Price.xlsx",
col_types = c("date", "numeric"))
Piggy <- Forecast_Bacon_Price
maybe <- Piggy %>%
mutate(Date = yearmonth(Month)) %>%
select(Date, Price) %>%
as_tsibble(index = Date)
pigplot <- autoplot(maybe) +
labs(y = "Price of Bacon (per Pound)", x = "Months", title = "Monthly Bacon Prices") +
geom_line(colour = "Black")## Plot variable not specified, automatically selected `.vars = Price`pigplot
pigstl <- maybe %>%
model(classical_decomposition(Price, type = "multiplicative")) %>%
components() %>%
autoplot() +
labs(title = "Multiplicative Decomp of Bacon")
pigstl## Warning: Removed 6 row(s) containing missing values (geom_path).
Accuracy
pigchoochoo <- maybe %>%
filter(year(Date) < 2016)
swineexam <- maybe %>%
filter(year(Date) > 2015)
pigparty <- pigchoochoo %>%
stretch_tsibble(.init = 10) %>%
model(
Holt = ETS(Price ~ error("A") + trend("A") + season("N")),
Component = ETS(Price ~ error("A") + trend("N") + season("N")),
Additive = ETS(Price ~ error("A") + trend("A") + season("A")),
Multiplicative = ETS(Price ~ error("M") + trend("A") + season("M")),
Damped = ETS(Price ~ error("A") + trend("Ad") + season("N"))) %>%
forecast(h = 4)## Warning: 8 errors (2 unique) encountered for Additive
## [3] A seasonal ETS model cannot be used for this data.
## [5] Not enough data to estimate this ETS model.## Warning: 8 errors (2 unique) encountered for Multiplicative
## [3] A seasonal ETS model cannot be used for this data.
## [5] Not enough data to estimate this ETS model.accuracy(pigparty, swineexam, list(rmse = RMSE, mae = MAE, mape = MAPE))## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing.
## 182 observations are missing between 2000 Nov and 2015 Dec## # A tibble: 5 × 5
## .model .type rmse mae mape
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 Additive Test 0.299 0.257 4.69
## 2 Component Test 0.297 0.260 4.74
## 3 Damped Test 0.396 0.352 6.41
## 4 Holt Test 0.349 0.317 5.77
## 5 Multiplicative Test 0.261 0.227 4.14Holt
Holtpig <- maybe %>%
model(AAN = ETS(Price ~ error("A") + trend("A") + season("N")))
augment(Holtpig) %>%
features(.resid, ljung_box, lag = 10, dof = 0)## # A tibble: 1 × 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 AAN 23.1 0.0105report(Holtpig)## Series: Price
## Model: ETS(A,A,N)
## Smoothing parameters:
## alpha = 0.9999
## beta = 0.0001103611
##
## Initial states:
## l[0] b[0]
## 2.499874 0.01890469
##
## sigma^2: 0.0161
##
## AIC AICc BIC
## 409.4910 409.7157 427.5384Holtpig %>%
forecast(h = 4) %>%
hilo()## # A tsibble: 4 x 6 [1M]
## # Key: .model [1]
## .model Date Price .mean `80%`
## <chr> <mth> <dist> <dbl> <hilo>
## 1 AAN 2022 Oct N(7.4, 0.016) 7.40 [7.239464, 7.564289]80
## 2 AAN 2022 Nov N(7.4, 0.032) 7.42 [7.191065, 7.650440]80
## 3 AAN 2022 Dec N(7.4, 0.048) 7.44 [7.158309, 7.720947]80
## 4 AAN 2023 Jan N(7.5, 0.064) 7.46 [7.133649, 7.783358]80
## # … with 1 more variable: `95%` <hilo>Component
Compopig <- maybe %>%
model(ETS(Price ~ error("A") + trend("N") + season("N")))
augment(Compopig) %>%
features(.resid, ljung_box, lag = 10, dof = 0)## # A tibble: 1 × 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 "ETS(Price ~ error(\"A\") + trend(\"N\") + season(\"N\"))" 22.6 0.0123report(Compopig)## Series: Price
## Model: ETS(A,N,N)
## Smoothing parameters:
## alpha = 0.9999
##
## Initial states:
## l[0]
## 2.750068
##
## sigma^2: 0.016
##
## AIC AICc BIC
## 407.0060 407.0952 417.8344Compopig %>%
forecast(h = 4) %>%
hilo()## # A tsibble: 4 x 6 [1M]
## # Key: .model [1]
## .model Date Price .mean `80%`
## <chr> <mth> <dist> <dbl> <hilo>
## 1 "ETS(Price ~ error(\"A\")… 2022 Oct N(7.4, 0.016) 7.38 [7.220737, 7.545261]80
## 2 "ETS(Price ~ error(\"A\")… 2022 Nov N(7.4, 0.032) 7.38 [7.153538, 7.612460]80
## 3 "ETS(Price ~ error(\"A\")… 2022 Dec N(7.4, 0.048) 7.38 [7.101972, 7.664026]80
## 4 "ETS(Price ~ error(\"A\")… 2023 Jan N(7.4, 0.064) 7.38 [7.058500, 7.707498]80
## # … with 1 more variable: `95%` <hilo>Additive
Pigplus <- maybe %>%
model(additive = ETS(Price ~ error("A") + trend("A") + season("A")))
augment(Pigplus) %>%
features(.resid, ljung_box, lag = 10, dof = 0)## # A tibble: 1 × 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 additive 15.5 0.116report(Pigplus)## Series: Price
## Model: ETS(A,A,A)
## Smoothing parameters:
## alpha = 0.9803652
## beta = 0.0001007609
## gamma = 0.013936
##
## Initial states:
## l[0] b[0] s[0] s[-1] s[-2] s[-3] s[-4]
## 3.051095 0.01350681 -0.09506619 -0.03456496 0.09526095 0.1464967 0.1563638
## s[-5] s[-6] s[-7] s[-8] s[-9] s[-10]
## 0.07218908 0.02523609 -0.03567011 -0.08330036 -0.0904351 -0.07813024
## s[-11]
## -0.07837957
##
## sigma^2: 0.0151
##
## AIC AICc BIC
## 404.9479 407.3479 466.3090Pigplus %>%
forecast(h = 4) %>%
hilo()## # A tsibble: 4 x 6 [1M]
## # Key: .model [1]
## .model Date Price .mean `80%`
## <chr> <mth> <dist> <dbl> <hilo>
## 1 additive 2022 Oct N(7.3, 0.015) 7.34 [7.186163, 7.501558]80
## 2 additive 2022 Nov N(7.2, 0.03) 7.23 [7.008111, 7.449813]80
## 3 additive 2022 Dec N(7.2, 0.044) 7.18 [6.911598, 7.450807]80
## 4 additive 2023 Jan N(7.2, 0.059) 7.20 [6.890578, 7.512198]80
## # … with 1 more variable: `95%` <hilo>Multiplicative
Pigmulti <- maybe %>%
model(multiplicative = ETS(Price ~ error("A") + trend("A") + season("A")))
augment(Pigmulti) %>%
features(.resid, ljung_box, lag = 10, dof = 0)## # A tibble: 1 × 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 multiplicative 15.5 0.116report(Pigmulti)## Series: Price
## Model: ETS(A,A,A)
## Smoothing parameters:
## alpha = 0.9803652
## beta = 0.0001007609
## gamma = 0.013936
##
## Initial states:
## l[0] b[0] s[0] s[-1] s[-2] s[-3] s[-4]
## 3.051095 0.01350681 -0.09506619 -0.03456496 0.09526095 0.1464967 0.1563638
## s[-5] s[-6] s[-7] s[-8] s[-9] s[-10]
## 0.07218908 0.02523609 -0.03567011 -0.08330036 -0.0904351 -0.07813024
## s[-11]
## -0.07837957
##
## sigma^2: 0.0151
##
## AIC AICc BIC
## 404.9479 407.3479 466.3090Pigmulti %>%
forecast(h = 4) %>%
hilo()## # A tsibble: 4 x 6 [1M]
## # Key: .model [1]
## .model Date Price .mean `80%`
## <chr> <mth> <dist> <dbl> <hilo>
## 1 multiplicative 2022 Oct N(7.3, 0.015) 7.34 [7.186163, 7.501558]80
## 2 multiplicative 2022 Nov N(7.2, 0.03) 7.23 [7.008111, 7.449813]80
## 3 multiplicative 2022 Dec N(7.2, 0.044) 7.18 [6.911598, 7.450807]80
## 4 multiplicative 2023 Jan N(7.2, 0.059) 7.20 [6.890578, 7.512198]80
## # … with 1 more variable: `95%` <hilo>Damped
Damnpig <- maybe %>%
model(`Damped Holt's method` = ETS(Price ~ error("A") + trend("Ad", phi = 0.9) + season("N")))
augment(Damnpig) %>%
features(.resid, ljung_box, lag = 10, dof = 0)## # A tibble: 1 × 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 Damped Holt's method 23.9 0.00793report(Damnpig)## Series: Price
## Model: ETS(A,Ad,N)
## Smoothing parameters:
## alpha = 0.9998997
## beta = 0.06998557
## phi = 0.9
##
## Initial states:
## l[0] b[0]
## 2.546801 0.08026293
##
## sigma^2: 0.0161
##
## AIC AICc BIC
## 409.7725 409.9972 427.8198Damnpig %>%
forecast(h = 4) %>%
hilo()## # A tsibble: 4 x 6 [1M]
## # Key: .model [1]
## .model Date Price .mean `80%`
## <chr> <mth> <dist> <dbl> <hilo>
## 1 Damped Holt's method 2022 Oct N(7.4, 0.016) 7.39 [7.231083, 7.556681]80
## 2 Damped Holt's method 2022 Nov N(7.4, 0.034) 7.40 [7.166093, 7.641258]80
## 3 Damped Holt's method 2022 Dec N(7.4, 0.055) 7.41 [7.113046, 7.711933]80
## 4 Damped Holt's method 2023 Jan N(7.4, 0.077) 7.42 [7.065480, 7.775365]80
## # … with 1 more variable: `95%` <hilo>Plot of All
pigpig <- maybe %>%
model(
Holt = ETS(Price ~ error("A") + trend("A") + season("N")),
Component = ETS(Price ~ error("A") + trend("N") + season("N")),
Additive = ETS(Price ~ error("A") + trend("A") + season("A")),
Multiplicative = ETS(Price ~ error("A") + trend("A") + season("A")),
`Damped Holt's method` = ETS(Price ~ error("A") + trend("Ad", phi = 0.9) + season("N"))) %>%
forecast(h = 4)
pigpig %>%
autoplot(maybe) +
geom_line(aes(y = .fitted), color = "red", data = augment(Damnpig)) +
geom_line(aes(y = .fitted), color = "blue", data = augment(Holtpig)) +
geom_line(aes(y = .fitted), color = "green", data = augment(Pigplus)) +
geom_line(aes(y = .fitted), color = "yellow", data = augment(Pigmulti)) +
labs(title = "Plot of all the methods")
Pigmulti %>%
forecast(h = 4) %>%
autoplot(maybe) +
geom_line(aes(y = .fitted), color = "red", data = augment(Pigmulti)) +
labs(y = "Bacon Price")
The Holt-Winter’s Multiplicative method seems to be my best option. The data isn’t lagged behind the original line, and I like that forecast shows the price of bacon falling in the future. The trend seems to say otherwise, but the large variation of prices suggest that bacon is somewhat due to fall soon. Using the ljung box test I can say that the residuals are statistically significant from zero.(Which makes me believe the forecast more) Although other models have a better fit, I believe this method is best mostly based on the AICc being much lower than the other methods. Finally this conclusion makes a lot of sense because the multiplicative method is used when the variation in the seasonal pattern, or the variation around the trend-cycle, appears to be proportional to the level of the time series. There are only two “shocks” in this model where this is untrue. (Early 2010 and Mid 2020)