Do exercises 8.1, 8.5, 8.6, 8.7, 8.8, 8.9 in Hyndman. Please submit both the link to your Rpubs and the .pdf file with your run code

library("fpp3")
## Registered S3 method overwritten by 'tsibble':
##   method               from 
##   as_tibble.grouped_df dplyr
## ── Attaching packages ──────────────────────────────────────────── fpp3 1.0.1 ──
## ✔ tibble      3.2.1     ✔ tsibble     1.1.6
## ✔ dplyr       1.1.4     ✔ tsibbledata 0.4.1
## ✔ tidyr       1.3.1     ✔ feasts      0.4.1
## ✔ lubridate   1.9.3     ✔ fable       0.4.1
## ✔ ggplot2     3.5.1
## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## ✖ lubridate::date()    masks base::date()
## ✖ dplyr::filter()      masks stats::filter()
## ✖ tsibble::intersect() masks base::intersect()
## ✖ tsibble::interval()  masks lubridate::interval()
## ✖ dplyr::lag()         masks stats::lag()
## ✖ tsibble::setdiff()   masks base::setdiff()
## ✖ tsibble::union()     masks base::union()
library("tidyverse")
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0     ✔ readr   2.1.5
## ✔ purrr   1.0.2     ✔ stringr 1.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter()     masks stats::filter()
## ✖ tsibble::interval() masks lubridate::interval()
## ✖ dplyr::lag()        masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library("ggplot2")

Question 8.1

Consider the the number of pigs slaughtered in Victoria, available in the aus_livestock dataset.

a.) Use the ETS() function to estimate the equivalent model for simple exponential smoothing. Find the optimal values of
α and ℓ0, and generate forecasts for the next four months.

data("aus_livestock")
aus_livestock
## # A tsibble: 29,364 x 4 [1M]
## # Key:       Animal, State [54]
##       Month Animal                     State                        Count
##       <mth> <fct>                      <fct>                        <dbl>
##  1 1976 Jul Bulls, bullocks and steers Australian Capital Territory  2300
##  2 1976 Aug Bulls, bullocks and steers Australian Capital Territory  2100
##  3 1976 Sep Bulls, bullocks and steers Australian Capital Territory  2100
##  4 1976 Oct Bulls, bullocks and steers Australian Capital Territory  1900
##  5 1976 Nov Bulls, bullocks and steers Australian Capital Territory  2100
##  6 1976 Dec Bulls, bullocks and steers Australian Capital Territory  1800
##  7 1977 Jan Bulls, bullocks and steers Australian Capital Territory  1800
##  8 1977 Feb Bulls, bullocks and steers Australian Capital Territory  1900
##  9 1977 Mar Bulls, bullocks and steers Australian Capital Territory  2700
## 10 1977 Apr Bulls, bullocks and steers Australian Capital Territory  2300
## # ℹ 29,354 more rows
help("aus_livestock")
pig_meat  = aus_livestock %>% 
  filter(State == "Victoria", Animal == "Pigs")
pig_meat |>
  autoplot(Count) +
  labs(y = "Pigs Slaughtered", title = "Pig Production in Victoria")

# Use the ETS function to find the optimal values of alpha and level values , # SES = Simple Exponential Smoothing
pig_fit <- pig_meat |>
  model(SES = ETS(Count ~ error("A") + trend("N") +  season("N")))
# ^ help us find the model
# use report to get the values from the model 
model_report = pig_fit %>% 
  report()
## Series: Count 
## Model: ETS(A,N,N) 
##   Smoothing parameters:
##     alpha = 0.3221247 
## 
##   Initial states:
##      l[0]
##  100646.6
## 
##   sigma^2:  87480760
## 
##      AIC     AICc      BIC 
## 13737.10 13737.14 13750.07

In the report we find the optimal valiues for alpha is 0.3221247 and l[0] is 100646.6

To find the next four months or four periods

pig_forecast = pig_fit %>% 
  forecast(  h = 4) # four months ahead
# we are looking for the next four months so i can begin on later data like 2016
pig_forecast %>%
  autoplot(pig_meat %>%  filter(Month >= yearmonth("2016 Jan"))) %>%
  labs(y = "Pigs Slaughtered", title = "Pig Production in Victoria")
## [[1]]

## 
## $y
## [1] "Pigs Slaughtered"
## 
## $title
## [1] "Pig Production in Victoria"
## 
## attr(,"class")
## [1] "labels"

b.) Compute a 95% prediction interval for the first forecast using ^ y ± 1.96s where s is the standard deviation of the residuals. Compare your interval with the interval produced by R.

# first need get the residuals to calculate the standard deviation
pig_residuals =  residuals(pig_fit) 
pig_sd = sd(pig_residuals$.resid,na.rm = TRUE)
num1_pig_forecast = pig_forecast$.mean[1] # get the first pig forecast, in the first row
# we can can find  the plus and minus of the 1.96 confidence interval , and multiplying it by the standard deviations
upper_bound_pig = pig_forecast$.mean[1] - (pig_sd * 1.96)
lower_bound_pig = pig_forecast$.mean[1] + (pig_sd * 1.96)
print(upper_bound_pig)
## [1] 76871.01
print(lower_bound_pig)
## [1] 113502.1

Now we have to compare prediction intervals

# The hilo() function converts the forecast distributions into intervals (Chapter 5.5)
pig_pred_interval = hilo(pig_forecast, level = 95) # confidence level in 95 percent

In the pig meat production in the state of Victoria, in Jan 2019, the production of pig within the percent confidence interval is 76871 to 113502 pig produced.

Question 8.5

Data set global_economy contains the annual Exports from many countries. Select one country to analyse.

data("global_economy")
global_economy
## # A tsibble: 15,150 x 9 [1Y]
## # Key:       Country [263]
##    Country     Code   Year         GDP Growth   CPI Imports Exports Population
##    <fct>       <fct> <dbl>       <dbl>  <dbl> <dbl>   <dbl>   <dbl>      <dbl>
##  1 Afghanistan AFG    1960  537777811.     NA    NA    7.02    4.13    8996351
##  2 Afghanistan AFG    1961  548888896.     NA    NA    8.10    4.45    9166764
##  3 Afghanistan AFG    1962  546666678.     NA    NA    9.35    4.88    9345868
##  4 Afghanistan AFG    1963  751111191.     NA    NA   16.9     9.17    9533954
##  5 Afghanistan AFG    1964  800000044.     NA    NA   18.1     8.89    9731361
##  6 Afghanistan AFG    1965 1006666638.     NA    NA   21.4    11.3     9938414
##  7 Afghanistan AFG    1966 1399999967.     NA    NA   18.6     8.57   10152331
##  8 Afghanistan AFG    1967 1673333418.     NA    NA   14.2     6.77   10372630
##  9 Afghanistan AFG    1968 1373333367.     NA    NA   15.2     8.90   10604346
## 10 Afghanistan AFG    1969 1408888922.     NA    NA   15.0    10.1    10854428
## # ℹ 15,140 more rows

I will be using the country of Ghana

ghana_export <- global_economy |>
  filter(Country == "Ghana")

a.) Plot the Exports series and discuss the main features of the data.

head(ghana_export)
## # A tsibble: 6 x 9 [1Y]
## # Key:       Country [1]
##   Country Code   Year         GDP Growth      CPI Imports Exports Population
##   <fct>   <fct> <dbl>       <dbl>  <dbl>    <dbl>   <dbl>   <dbl>      <dbl>
## 1 Ghana   GHA    1960 1217230038.  NA    NA          35.4    28.2    6652287
## 2 Ghana   GHA    1961 1302674264.   3.43 NA          36.5    26.1    6866539
## 3 Ghana   GHA    1962 1382515590.   4.11 NA          28.5    24.2    7085464
## 4 Ghana   GHA    1963 1540797517.   4.41 NA          27.5    21.2    7303432
## 5 Ghana   GHA    1964 1731296119.   2.21  0.00131    23.9    19.9    7513289
## 6 Ghana   GHA    1965 2053462872.   1.37  0.00166    26.7    17.1    7710549
ghana_export %>%
  autoplot(Exports)  +
  labs(title = 'Ghana Export Trends: From 1960 to 2017', subtitle = "Annual Exoport Value Over Time")

This graph shows the Ghana exports form 1960 to 2000. The exports started around 28 units and steadily dropped throughout the years until the 1980s. From the mid 1980s to the year 2000, exports grew steadily, where it hit its peak of 50 units, with some fluctuations throughout the growth. The Ghana exports relies most on commodities like gold, oil and cocoa. These commodities are heavily impact by global market polices like tariffs and interest rates. After the year 2000, there was sharp decline, with the 2008 recession impact the economy. However, there there was been slow recovery on the export of commodities.

b.) Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.

# using the ANN model to forecast the series
gh_fit = ghana_export %>% 
  model(SES = ETS(Exports ~ error("A") + trend("N") + season("N")))
gh_ex_report = gh_fit %>% 
  report()
## Series: Exports 
## Model: ETS(A,N,N) 
##   Smoothing parameters:
##     alpha = 0.9755569 
## 
##   Initial states:
##      l[0]
##  28.14357
## 
##   sigma^2:  20.1376
## 
##      AIC     AICc      BIC 
## 413.6206 414.0650 419.8019
gh_forecast = gh_fit %>% 
  forecast( h = 8) # 8 years forecasted
gh_forecast %>% 
  autoplot(ghana_export) +  labs(title = 'Ghana Export Forecast For 8 Years', subtitle = "Annual Exoport Value Over Time")

b.) Compute the RMSE values for the training data.

gh_accuracy =accuracy(gh_fit)  # We use accuracy function RSME values
print(gh_accuracy) # RSME to help find how far the predicted values are from the actual vlaues
## # 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 Ghana   SES    Training 0.124  4.41  3.10 -2.10  16.6 0.986 0.991 0.000306
gh_accuracy %>% 
  pull(RMSE) # pull out the RMSE value
## [1] 4.409445

d.) Compare the results to those from an ETS(A,A,N) model. (Remember that the trended model is using one more parameter than the simpler model.) Discuss the merits of the two forecasting methods for this data set.

gh_aan_model = ghana_export %>% 
    model( Additive = ETS(Exports ~ error("A") + trend("A") + season("N")))
accuracy(gh_aan_model) # The ANN model is the additive used for constant variance, also called the Holts Linear Trend Model
## # 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 Ghana   Additive Training 0.0154  4.41  3.12 -2.80  16.8 0.993 0.991 -0.000769
gh_aan_accuracy = accuracy(gh_aan_model)

Comparing them both

print(gh_accuracy)
## # 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 Ghana   SES    Training 0.124  4.41  3.10 -2.10  16.6 0.986 0.991 0.000306
print(gh_aan_accuracy)
## # 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 Ghana   Additive Training 0.0154  4.41  3.12 -2.80  16.8 0.993 0.991 -0.000769

The Simple Exponential Smoothing Model has a lower RSME (4.4094) than our Additive Model of 4.4106 indicating it will have more accurate prediction of the export trend for the next 8 years. The ANN model observes trend for constant increase and decrease, however our original data show huge magnitude of fluctuations throughout the data, an AAN will not be useful for forecasting. Our SES model does perform better than the ANN model , because of our weak and strong our trend develop on sharp and low levels.

e.) Compare the forecasts from both methods. Which do you think is best?

gh_forecast_sse = gh_fit %>% 
  forecast(h = 8) # Using 8 years in the future for forecast
gh_forecast_aan = gh_aan_model %>%  # forecast for addtive model 
  forecast(h = 8)
autoplot(ghana_export, Exports) +
  autolayer(gh_forecast_sse, series = "Simple Exponential Smoothing   (ANN Model)", color = "blue") +
  autolayer(gh_forecast_aan, series = "Additive  Model (AAN)", color = "red") +
  labs(
    title = "Ghana Export Forecast in 8 Years)",
    subtitle = " SES vs Holts Additive Models",
    y = "Exports"
  ) 
## Warning in ggdist::geom_lineribbon(without(intvl_mapping, "colour_ramp"), :
## Ignoring unknown parameters: `series`
## Warning in geom_line(mapping = without(mapping, "shape"), data =
## unpack_data(object[single_row[["FALSE"]], : Ignoring unknown parameters:
## `series`
## Warning in ggdist::geom_lineribbon(without(intvl_mapping, "colour_ramp"), :
## Ignoring unknown parameters: `series`
## Warning in geom_line(mapping = without(mapping, "shape"), data =
## unpack_data(object[single_row[["FALSE"]], : Ignoring unknown parameters:
## `series`
## Scale for fill_ramp is already present.
## Adding another scale for fill_ramp, which will replace the existing scale.

f.) Calculate a 95% prediction interval for the first forecast for each model, using the RMSE values and assuming normal errors. Compare your intervals with those produced using R.

num1_sse_forecast = gh_forecast %>%  
  slice(1) %>% # get the first row in the SES Model
 pull(.mean) # get the mean value 
nuum1_add_forecast = gh_forecast_aan %>% 
  slice(1) %>%  # get the first row for additive model 
  pull(.mean)
rmse_sse = gh_accuracy  %>% 
  pull(RMSE)
rmse_add = gh_aan_accuracy %>% 
  pull(RMSE)
# Simple Smoothing Prediction Interval in the 95 percent confidence interval 
lower_sse <- num1_sse_forecast - 1.96 * rmse_sse
upper_sse <- num1_sse_forecast+ 1.96 * rmse_sse
# Holts Linear Prediction Interval 
lower_add <- nuum1_add_forecast - 1.96 * rmse_add
upper_add <- nuum1_add_forecast + 1.96 * rmse_add
# we can use the "cat" to show the results
cat("Simple Exponential Smoothing ETS Method 95% Confidence Interval: [", lower_sse, ",", upper_sse, "]\n")
## Simple Exponential Smoothing ETS Method 95% Confidence Interval: [ 26.53294 , 43.81796 ]
cat("Holt's Additive Trend in 95%  Confidence  Interval: [", lower_add, ",", upper_add, "]\n")
## Holt's Additive Trend in 95%  Confidence  Interval: [ 26.65874 , 43.94842 ]

Question 8.6 Forecast the Chinese GDP from the global_economy data set using an ETS model. Experiment with the various options in the ETS() function to see how much the forecasts change with damped trend, or with a Box-Cox transformation. Try to develop an intuition of what each is doing to the forecasts.

china_gdp = global_economy %>% 
  filter(Country == 'China')
china_gdp %>% autoplot(GDP) +
  labs(title = 'China Gross Domestic Product')

Experiment with the various options in the ETS() function to see how much the forecasts change with damped trend

china_sse_fit = china_gdp %>% 
  model(
    SSE = ETS(GDP ~ error("A") + trend("N") + season("N")))
china_holt_fit = china_gdp %>% 
  model(Holt = ETS(GDP ~ error("A") + trend("A") + season("N")))
china_damped_fit = china_gdp %>% 
  model(
    Damped = ETS(GDP ~ error("A") + trend("Ad") + season("N"))
  )
# Forecasting for the next 15 years 
china_forecast_sse = china_sse_fit %>% 
  forecast(h = 15) 
china_forecast_holt = china_holt_fit %>% 
  forecast(h = 15)
china_forecast_damp = china_damped_fit %>% 
  forecast(h = 15)
china_forecast_sse %>% 
  autoplot(china_gdp) + autolayer(china_forecast_holt, series = "Holt", colour = "yellow") + autolayer(china_forecast_damp, series = "Damped", color = "red")  +
  
  theme_update() + labs(title = "China GDP Forecast For 8 Years", subtitle = "Comparison of ETS Model: Simple Exponential Smoothing, Holt's Linear Trend, and Damped Models ") 
## Warning in ggdist::geom_lineribbon(without(intvl_mapping, "colour_ramp"), :
## Ignoring unknown parameters: `series`
## Warning in geom_line(mapping = without(mapping, "shape"), data =
## unpack_data(object[single_row[["FALSE"]], : Ignoring unknown parameters:
## `series`
## Scale for fill_ramp is already present.
## Adding another scale for fill_ramp, which will replace the existing scale.
## Warning in ggdist::geom_lineribbon(without(intvl_mapping, "colour_ramp"), : Ignoring unknown parameters: `series`
## Ignoring unknown parameters: `series`
## Scale for fill_ramp is already present.
## Adding another scale for fill_ramp, which will replace the existing scale.

Question 8.7

Find an ETS model for the Gas data from aus_production and forecast the next few years. Why is multiplicative seasonality necessary here? Experiment with making the trend damped. Does it improve the forecasts?

data("aus_production")
aus_production
## # A tsibble: 218 x 7 [1Q]
##    Quarter  Beer Tobacco Bricks Cement Electricity   Gas
##      <qtr> <dbl>   <dbl>  <dbl>  <dbl>       <dbl> <dbl>
##  1 1956 Q1   284    5225    189    465        3923     5
##  2 1956 Q2   213    5178    204    532        4436     6
##  3 1956 Q3   227    5297    208    561        4806     7
##  4 1956 Q4   308    5681    197    570        4418     6
##  5 1957 Q1   262    5577    187    529        4339     5
##  6 1957 Q2   228    5651    214    604        4811     7
##  7 1957 Q3   236    5317    227    603        5259     7
##  8 1957 Q4   320    6152    222    582        4735     6
##  9 1958 Q1   272    5758    199    554        4608     5
## 10 1958 Q2   233    5641    229    620        5196     7
## # ℹ 208 more rows
gas_prod = aus_production %>% 
  select(Quarter,Gas)
gas_prod %>% 
  autoplot()
## Plot variable not specified, automatically selected `.vars = Gas`

gas_fit = gas_prod %>% 
  model(ETS(Gas))
report(gas_fit)
## Series: Gas 
## Model: ETS(M,A,M) 
##   Smoothing parameters:
##     alpha = 0.6528545 
##     beta  = 0.1441675 
##     gamma = 0.09784922 
## 
##   Initial states:
##      l[0]       b[0]      s[0]    s[-1]    s[-2]     s[-3]
##  5.945592 0.07062881 0.9309236 1.177883 1.074851 0.8163427
## 
##   sigma^2:  0.0032
## 
##      AIC     AICc      BIC 
## 1680.929 1681.794 1711.389
gas_fit_forecast = gas_fit %>% 
  forecast(h = 12) 
gas_fit_forecast %>% 
  autoplot(gas_prod  %>%  filter(Quarter >= yearquarter("2005 Q1")))

Now with multiplicative seasonality ^ same model from before ( Rewrote code here)

gas_multi_fit = gas_prod %>% 
  model( Mulit = ETS(Gas ~ error("M") + trend("A") + season("M")))
gas_multi_fit_forecast = gas_multi_fit %>% 
  forecast(h = 12) 
gas_multi_fit_forecast %>% 
  autoplot(gas_prod  %>%  filter(Quarter >= yearquarter("2005  Q1"))) + labs(title = "Multiplicative Seasonality Model for Gas Production Forecasting ")

Experimenting with the damped additive model

gas_damped_fit = gas_prod %>% 
  model( Damped = ETS(Gas ~ error("A") + trend("Ad") + season("M"))) # Ad for  damped, trend can gradually decrease or become stable
gas_forecast_damped <- gas_damped_fit  %>% 
  forecast(h = 12)
gas_forecast_damped %>% 
  autoplot(gas_prod  %>%  filter(Quarter >= yearquarter("2005  Q1"))) + labs(title = "Damped Additive Model for Gas Production Forecasting")

Multiplicative Seasonality is necessary here due the seasonal peak seen through out the quarters. The demand for gas is most dependent on the month the consumers are using their gas powered product. The seasonal fluctuations are observed to increased proportionally in the series, and multiplicative model will capture trends increasing in a proportional level.However, using the damped additive with a multiplicative seasonal component shows that we expect the demand for gas to stay stable overtime. This approach to observing the data, is not useful as in the future demand for gas been influenced by factors such as rise in population, construction of new houses, and technoligcal advances..

Question 8.8

Recall your retail time series data (from Exercise 7 in Section 2.10).

set.seed(12345678)
myseries <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`,1))
myseries %>% 
  autoplot()
## Plot variable not specified, automatically selected `.vars = Turnover`

a.) Why is multiplicative seasonality necessary for this series?

We use multiplicative seasonality due to the seasonal turnover fluctuation increasing throughout the years. There are peaks in which people are hired for seasonal positions which a retial locaiton expects an high increase of demand. This months are summer and the winter, combined with popular holidays.

b.) Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.

retail_multi_fit = myseries %>% 
    model( Holt_Winters = ETS(Turnover ~ error("M") + trend("A") + season("M")),
        Damped_Winters =  ETS(Turnover ~ error("M") + trend("Ad") + season("M"))
                    )
retail_holt_forecast = retail_multi_fit %>% 
  forecast(h = 12)
retail_holt_forecast %>% 
  autoplot(myseries) + labs(title = " Turnover in Retail", subtitle = "Holt-Winters vs  Damped Winters Models" )

c.) Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?

retail_accuracy = accuracy(retail_multi_fit)
view(retail_accuracy)

Comparing our RMSE values in the one step forecast with the Holt Winters and Damped Winters, we get the value 0.6130408 for Holts Winter and 0.6155067 for Damped Winters model. The Holts Winter model has lower RMSE, indicating a better fit for accuary in the forecasting data.

d.) Check that the residuals from the best method look like white noise.

retail_multi_fit %>% 
  select(Holt_Winters) %>% 
   gg_tsresiduals()

e.) Now find the test set RMSE, while training the model to the end of 2010. Can you beat the seasonal naïve approach from Exercise 7 in Section 5.11?

myseries2010 = myseries %>% 
  filter(year(Month) < 2011)
myseries2010_fit = myseries2010 %>% 
    model( Holt_Winters = ETS(Turnover ~ error("M") + trend("A") + season("M")),
        Damped_Winters =  ETS(Turnover ~ error("M") + trend("Ad") + season("M")),
        Seasonal_Naive = SNAIVE(Turnover) # seasonal naive method from chapter 5
                    )
myseries2010_forecast = myseries2010_fit %>%
  forecast(new_data = anti_join(myseries, myseries2010)) # using anti join to filer out rows that been in myseries already, make the data much cleaner
## Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover)`
myseries2010_forecast %>%
  autoplot(myseries %>%  filter(Month >= yearmonth("2010 Jan")), level = NULL) +
  labs(y = "Turnover", title = "2010 and Beyond Retial  Turnover Forecasting ", subtitle = "Holt_Winters  vs  Damped_Winters vs  Seasonal_Naive  ") +
  guides(color = guide_legend(title = "Model"))

Question 8.9 For the same retail data, try an STL decomposition applied to the Box-Cox transformed series, followed by ETS on the seasonally adjusted data. How does that compare with your best previous forecasts on the test set?

series_lambda <- myseries2010 |>
  features(Turnover, features = guerrero) |> # function to find lambda
  pull(lambda_guerrero)
series_boxcox = myseries2010 %>% 
  mutate(
    Boxcox = box_cox(Turnover,series_lambda)
  )  #  make row for the boxcox of each turnover, box used for normal distribtion  of data
series_boxcox_fit <- series_boxcox |> 
  model(
    ETS = ETS(Boxcox), # find the best compoenents
    Seasonal_Component = STL(Boxcox ~ season(window = "periodic")) # for decomposing the compoents
  )
# Graph looked funky b4, fixed now adjusted with element text(horizonal adjustment)
series_boxcox_fit |> 
  components() |> 
  autoplot() +
  labs(
    title = "Decomposition of Box-Cox Transformed Series",
    subtitle = "Trend, Seasonal, and Remainder Components of Boxcox ",
    y = "Turnover Boxcox",
    x = "Month"
  ) +
  theme_minimal() +  # Clean theme
  theme(
    plot.title = element_text(hjust = 0.5), # bring the title to middle of the graph
    plot.subtitle = element_text(hjust = 0.5),  
    legend.position = "bottom"  #  move legend from the side to the bottom to looker better
  ) +
  scale_color_brewer(palette = "Set2") #ggplot has lots fo color palette so this looks
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).

box_cox_hw = myseries2010 %>% 
  model(
    MAN_model = ETS(Turnover ~ error("M") + trend("A") + season("N")))
box_cox_hw_forecast = box_cox_hw %>% 
  forecast(h = 12) 
box_cox_hw_forecast %>% 
  autoplot(myseries2010) + labs(title = "Retial Turnover Boxcox Forecast for 12 Months")