library(readxl)
library(tidyquant)
library(fpp3)
library(moments)
library(tsibble)
library(tsibbledata)
library(ggplot2)
library(dplyr)
library(seasonal)
library(seasonalview)
library(fabletools)
library(fable)
library(ggfortify)
library(quantmod)
library(GGally)
Bacon <- read_excel("C:/Users/harle/OneDrive/Desktop/Class Files/Fall 2022/Forecasting/Bacon.xlsx")

bacon_tsibble <- Bacon %>%
  mutate(Month = yearmonth(Month)) %>%
  as_tsibble(index = Month)
autoplot(bacon_tsibble) +
  labs(title = "Monthly Price of Bacon, 1/2000-9/2022")
## Plot variable not specified, automatically selected `.vars = Price`

#There is a very clear upward trend but no cyclicity or evident seasonality.

exponential_smooth <- bacon_tsibble %>%
  model(ETS(Price ~ error("A") + trend("N") + season("N")))
  
exponential_smooth %>%
  forecast(h = 4)%>%
  autoplot(bacon_tsibble) +
  geom_line(aes(y = .fitted), col="#D55E00",
            data = augment(exponential_smooth)) +
  labs(title = "Four Month Forecast of Bacon Prices with Simple Smooth")

#This method should not be used since the data has an obvious upward trend. The simple method should only be used when there is no trend or seasonality evident in the initial data.

bacon_tsibble %>%
  model(
    `Holt's method` = ETS(Price ~ error("A") +  trend("A") + season("N"))) %>%
  forecast(h = 4) %>%
  autoplot(bacon_tsibble) +
  labs(title = "Four Month Forecast of Bacon Prices with Holt's Linear Trend")

#This method should yield a relatively accurate forecast, as it is the simple method plus a trend component. I do not believe a damped method would be more accurate as overforecasting should not be an issue since the forecasted period is not very far into the future. 

bacon_tsibble %>%
  model(
    `Damped Holt's method` = ETS(Price ~ error("A") +
                       trend("Ad", phi = 0.9) + season("N"))) %>%
  forecast(h = 4) %>%
  autoplot(bacon_tsibble) +
  labs(title = "Four Month Forecast of Bacon Prices with Holt's Damped Linear Trend")

#To test my theory about the damped method not being more accurate I produced a forecast with Holt's Damped Linear Trend. This forecast and the regular Holt's Linear Trend forecast are nearly identical. I suspect this is because the requested forecast is only for four months into the future. This method and the HLT method should have very close RMSE and MAPE scores.

bacon_tsibble %>%
  model(
    additive = ETS(Price ~ error("A") + trend("A") + season("A"))) %>%
  forecast(h = 4) %>%
  autoplot(bacon_tsibble) +
  labs(title = "Four Month Forecast of Bacon Prices with Holt-Winter's Additive Method")

#This method adds a seasonal component to Holt's Linear Trend. The additive method is preferable when there is constant seasonality. I did not note any seasonality when I first plotted, so I do not expect any additional accuracy from this method.

bacon_tsibble %>%
  model(
    multiplicative = ETS(Price ~ error("M") + trend("A") + season("M"))) %>%
  forecast(h = 4) %>%
  autoplot(bacon_tsibble) +
  labs(title = "Four Month Forecast of Bacon Prices with Holt-Winter's Multiplicative Method")

#Checking out multiplicative just to see. One should use multiplicative when seasonality is not constant. I did not note any seasonality in my first plot so like with HW's additive I do not expect any additional degree of accuracy.

bacon_tsibble %>%
  stretch_tsibble(.init = 10) %>%
  model(
    SES = ETS(Price ~ error("A") + trend("N") + season("N")),
    Holt = ETS(Price ~ error("A") + trend("A") + season("N")),
    Damped = ETS(Price ~ error("A") + trend("Ad") + season("N")),
    multiplicative = ETS(Price ~ error("M") + trend("A") + season("M")),
    additive = ETS(Price ~ error("A") + trend("A") + season("A"))
  ) %>%
  forecast(h = 1) %>%
  accuracy(bacon_tsibble)
## 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.
## 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: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing. 
## 1 observation is missing at 2022 Oct
## # A tibble: 5 × 10
##   .model         .type       ME  RMSE    MAE     MPE  MAPE  MASE RMSSE  ACF1
##   <chr>          <chr>    <dbl> <dbl>  <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 additive       Test   0.00148 0.135 0.0997 -0.0440  2.16 0.289 0.272 0.320
## 2 Damped         Test   0.0132  0.129 0.0944  0.248   2.04 0.274 0.260 0.242
## 3 Holt           Test   0.00176 0.128 0.0949 -0.0315  2.06 0.275 0.259 0.255
## 4 multiplicative Test  -0.00120 0.142 0.107  -0.118   2.33 0.311 0.286 0.376
## 5 SES            Test   0.0166  0.128 0.0941  0.304   2.02 0.273 0.258 0.248
#The SES has the lowest RMSE, MAE, AND MAPE scores so I will use it for my forecast. 
bacon_forecast <- bacon_tsibble %>%
  model(ETS(Price ~ error("A") + trend("N") + season("N"))) %>%
  forecast(h=4)

bacon_forecast %>%
  hilo(80)
## # A tsibble: 4 x 5 [1M]
## # Key:       .model [1]
##   .model                        Month         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
#My point forecast is $7.38 with an 80% prediction interval of [7.058500, 7.707498].