library(dplyr)
library(distributional)
library(fable)
library(feasts)
library(forecast)
library(fpp3)
library(ggfortify)
library(httr)
library(lubridate)
library(readr)
library(readxl)
library(seasonal)
library(stats)
library(tsibble)
library(tsibbledata)
library(tidyr)

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.

f <- aus_livestock %>%
  filter(State == "Victoria",
         Animal == "Pigs") %>%
  model(ETS(Count ~ error("A") + trend("N") + season("N")))
f_forecast <- f %>%
  forecast(h = 4)
f_forecast %>%
  autoplot(aus_livestock)

report(f)
## 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
f_forecast
## # A fable: 4 x 6 [1M]
## # Key:     Animal, State, .model [1]
##   Animal State    .model                          Month             Count  .mean
##   <fct>  <fct>    <chr>                           <mth>            <dist>  <dbl>
## 1 Pigs   Victoria "ETS(Count ~ error(\"A\") +… 2019 Jan N(95187, 8.7e+07) 95187.
## 2 Pigs   Victoria "ETS(Count ~ error(\"A\") +… 2019 Feb N(95187, 9.7e+07) 95187.
## 3 Pigs   Victoria "ETS(Count ~ error(\"A\") +… 2019 Mar N(95187, 1.1e+08) 95187.
## 4 Pigs   Victoria "ETS(Count ~ error(\"A\") +… 2019 Apr N(95187, 1.1e+08) 95187.
  1. 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.
sd_residuals <- residuals(f)$.resid %>% sd()
hat_y <- f_forecast$.mean[1]
f_forecast %>% hilo(95) %>% pull('95%') %>% head(1)
## <hilo[1]>
## [1] [76854.79, 113518.3]95

8.5 Data set global_economy contains the annual Exports from many countries. Select one country to analyse.
a. Plot the Exports series and discuss the main features of the data.
We can see large increase during beginning of the 1980s and then it decreases in the 1990’s. There is no seasonality.

global_economy %>%
  filter(Country == "Turkey") %>%
  autoplot(Exports) +
  labs(y="% GDP")

  1. Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
turkeyFit <- global_economy %>%
  filter(Country == "Turkey") %>%
  model(ETS(Exports ~ error("A") + trend("N") + season("N"))) 

forecastTurkey <- turkeyFit %>%
  forecast(h = 5)

forecastTurkey %>%
  autoplot(global_economy) +
  labs(y="% GDP")

  1. Compute the RMSE values for the training data.
accuracy(turkeyFit)$RMSE
## [1] 2.183255
  1. 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.
    ETS(A,A,N) may have more optimal forecasting since its RMSE is lower.
turkeyFit2 <- global_economy %>%
  filter(Country == "Turkey") %>%
  model(ETS(Exports ~ error("A") + trend("A") + season("N")))

forecastTurkey2 <- turkeyFit2 %>%
  forecast(h = 4)

forecastTurkey2 %>%
  autoplot(global_economy) +
  labs(y="% GDP")

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

ETS(A,A,N) would be best because it displays lower RMSE and shows increasing trend in the data. The other method suggests does not show much change in the exports.

  1. 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.
    My code was not working properly to fully answer this question but I anticipate using RMSE as s, that the interval would be larger.
x <- residuals(turkeyFit)$.resid %>% sd()
y <- forecastTurkey$.mean[1]

forecastTurkey %>% hilo(95) %>% pull('95%') %>% head(1)
## <hilo[1]>
## [1] [20.14458, 28.85427]95
x2 <- residuals(turkeyFit2)$.resid %>% sd()
y2 <- forecastTurkey2$.mean[1]

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.
[Hint: use a relatively large value of h when forecasting, so you can clearly see the differences between the various options when plotting the forecasts.]
There is very little seasonality in China GDP.

l <- global_economy %>%
  filter(Country == "China") %>%
  features(GDP, features = guerrero) %>%
  pull(lambda_guerrero)

chinaFit <- global_economy %>%
  filter(Country == "China") %>%
  model(`Simple` = ETS(GDP ~ error("A") + trend("N") + season("N")),
        `Holt's method` = ETS(GDP ~ error("A") + trend("A") + season("N")),
        `Damped Holt's method` = ETS(GDP ~ error("A") + trend("Ad", phi = 0.8) + season("N")),
        `Box-Cox` = ETS(box_cox(GDP,l) ~ error("A") + trend("A") + season("N")),
        `Box-Cox Damped` = ETS(box_cox(GDP,l) ~ error("A") + trend("Ad", phi = 0.8) + season("N")),
        `Log` = ETS(log(GDP) ~ error("A") + trend("A") + season("N")),
        `Log Damped` = ETS(log(GDP) ~ error("A") + trend("Ad", phi = 0.8) + season("N"))
        ) 

forecastChina <- chinaFit %>%
  forecast(h = 14)

forecastChina %>%
  autoplot(global_economy, level = NULL)

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?
It is needed because the seasonal pattern’s variation seems proportional to the time series level. The amplitude increases with the trend increasing.

gas <- aus_production %>%
  model(additive = ETS(Gas ~ error("A") + trend("A") + season("A")),
        multiplicative = ETS(Gas ~ error("M") + trend("A") + season("M")),
        `damped multiplicative` = ETS(Gas ~ error("M") + trend("Ad", phi = 0.9) + season("M"))) 
aus_production %>%
   model(additive = ETS(Gas ~ error("A") + trend("A") + season("A")),
        multiplicative = ETS(Gas ~ error("M") + trend("A") + season("M"))) %>%
   forecast(h=19) %>%
  autoplot(aus_production, level = NULL)

aus_production %>%
  model(multiplicative = ETS(Gas ~ error("M") + trend("A") + season("M")),
        `damped multiplicative` = ETS(Gas ~ error("M") + trend("Ad", phi = 0.9) + season("M"))
        ) %>%
  forecast(h=19) %>%
  autoplot(aus_production, level= NULL)

report(gas)
## Warning in report.mdl_df(gas): Model reporting is only supported for individual
## models, so a glance will be shown. To see the report for a specific model, use
## `select()` and `filter()` to identify a single model.
## # A tibble: 3 × 9
##   .model                  sigma2 log_lik   AIC  AICc   BIC   MSE  AMSE    MAE
##   <chr>                    <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 additive              23.6       -927. 1872. 1873. 1903.  22.7  29.7 3.35  
## 2 multiplicative         0.00324   -831. 1681. 1682. 1711.  21.1  32.2 0.0413
## 3 damped multiplicative  0.00340   -835. 1688. 1689. 1719.  21.0  32.4 0.0424

8.8 Recall your retail time series data (from Exercise 8 in Section 2.10).
a. Why is multiplicative seasonality necessary for this series?
It is needed because the seasonal pattern’s variation seems proportional to the time series level. The amplitude increases with the trend increasing.
b. Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.

set.seed(1234)
series <- aus_retail %>%
  filter(`Series ID` == sample(aus_retail$`Series ID`,1)) 
curr <- series %>%
  model(multiplicative = ETS(Turnover ~ error("M") + trend("A") + season("M")),
        `damped multiplicative` = ETS(Turnover ~ error("M") + trend("Ad") + season("M"))) 
curr %>%
  forecast(h=36) %>%
  autoplot(series, level = NULL)

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

accuracy(curr) %>%
  select(.model, RMSE)
## # A tibble: 2 × 2
##   .model                 RMSE
##   <chr>                 <dbl>
## 1 multiplicative         1.34
## 2 damped multiplicative  1.36
  1. Check that the residuals from the best method look like white noise.
    The residuals like white noise. We conclude that we cannot distinguish the residuals from white noise series.
series %>%
  model(multiplicative = ETS(Turnover ~ error("M") + trend("A") + season("M"))) %>%
  gg_tsresiduals()

series %>%
  model(multiplicative = ETS(Turnover ~ error("M") + trend("A") + season("M"))) %>%
  augment() %>% 
  features(.innov, box_pierce, lag = 25, dof = 0)
## # A tibble: 1 × 5
##   State    Industry                                     .model bp_stat bp_pvalue
##   <chr>    <chr>                                        <chr>    <dbl>     <dbl>
## 1 Tasmania Cafes, restaurants and takeaway food servic… multi…    27.3     0.343
  1. 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?
    RMSE is lower compared to seasonal naive approach so the multiplicative method forecasts better.
trainSeries <- series %>%
  filter(year(Month) < 2011)
ft <- trainSeries %>%
  model(multi = ETS(Turnover ~ error("M") + trend("A") + season("M")),
        snaive = SNAIVE(Turnover))
forecast <- ft %>%
  forecast(new_data = anti_join(series, trainSeries))
## Joining, by = c("State", "Industry", "Series ID", "Month", "Turnover")
forecast %>% autoplot(series, level = NULL)

accuracy(ft) %>%
  select(.type, .model, RMSE)
## # A tibble: 2 × 3
##   .type    .model  RMSE
##   <chr>    <chr>  <dbl>
## 1 Training multi   1.18
## 2 Training snaive  2.90
forecast %>% accuracy(series)  %>%
  select(.type, .model, RMSE)
## # A tibble: 2 × 3
##   .type .model  RMSE
##   <chr> <chr>  <dbl>
## 1 Test  multi   3.99
## 2 Test  snaive  9.13

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?
My code was not working properly but I anticipate that the RMSE would be better for the test data.

#z <- forecast(stlf(trainSeries, lambda = BoxCox.lambda(trainSeries), h = 15), PI=TRUE)
#autoplot(z)