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

pigs <- aus_livestock|>
  filter(State == 'Victoria', Animal == 'Pigs')
autoplot(pigs, Count)+ 
  labs(y = "Count"
         , x = "Date"
         , title = "Pigs Slaughtered, Victoria")

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.

fit <- pigs|>
  model(ETS(Count ~ error("A")+trend("N")+season("N")))
report(fit)
## 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

According to the report of fit we can see that the optimal values for α and \(ℓ_0\) are .3221247 and 100646.6, respectively. Lets generate the next four forecast using the model.

fc <- fit |>
  forecast(h = 4)

The forecast are generated and let’s plot them on the real time series

pigs_new<- pigs %>%
  filter(yearmonth(Month) > yearmonth("2015 Dec")) %>%
  index_by(Date = as_date(Month))
fc|>
  autoplot(pigs_new)

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.

Let’s look at the intervals produced by R:

int <- hilo(fc)
print(int$`95%`[1])
## <hilo[1]>
## [1] [76854.79, 113518.3]95

Now let’s calculate manually for the first forecast

resids <- augment(fit)
s <- sd(resids$.resid) 
print(paste("Standard Deviation: ",s))
## [1] "Standard Deviation:  9344.66579258967"

We got the standard deviation of residuals now lets calculate the confidence interval

upper <- fc$.mean[1]+1.96*s
lower <- fc$.mean[1]-1.96*s

print(paste("Upper limit: ", upper))
## [1] "Upper limit:  113502.102384467"
print(paste("Lower limit: ", lower))
## [1] "Lower limit:  76871.0124775157"

As we can see that the values are really close to each other.

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

us_econ <- global_economy|>
  filter(Country == "United States", Year <= '2016')

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

us_econ|>
  autoplot(Exports)+labs(y= "Exports (% of GDP)", title = "Exports of United States in Percent of GDP")

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

Here is the ETS(A,N,N) model:

fit_us <- us_econ|>
  model(ETS(Exports ~ error("A")+trend("N")+season("N")))

Our model is ready let’s forecast the and in this we will forecast the next four years.

us_fc <- fit_us |>
  forecast(h = 4)

Now we have our forecast and let’s plot the forecast.

us_fc|>
  autoplot(us_econ)

c. Compute the RMSE values for the training data.

accuracy(fit_us)
## # 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 United States "ETS(Expo… Trai… 0.125 0.627 0.465  1.37  5.10 0.990 0.992 0.239

The RMSE value is .6270.

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.

mix_us <- us_econ %>%
  model(
    ANN= ETS(Exports ~error("A")+trend("N")+season("N"))
    ,AAN= ETS(Exports ~error("A")+trend("A")+season("N"))
    
    )

accuracy(mix_us) 
## # 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 United States ANN    Train… 0.125  0.627 0.465  1.37    5.10 0.990 0.992 0.239
## 2 United States AAN    Train… 0.0108 0.615 0.466 -0.0570  5.19 0.992 0.973 0.238

AS we can see that the RMSE for AAN model is lower than ANN.

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

mix_us %>%
  forecast(h = 5) %>%
  autoplot(us_econ, level=NULL) +
  labs(y = "Count", title = "Mexico exports with Forecast"
       , x = "Date")

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.

mix_us %>%
  forecast(h=1) %>%
  mutate(interval = hilo(Exports, 95))
## # A fable: 2 x 6 [1Y]
## # Key:     Country, .model [2]
##   Country       .model  Year     Exports .mean               interval
##   <fct>         <chr>  <dbl>      <dist> <dbl>                 <hilo>
## 1 United States ANN     2017 N(12, 0.41)  11.9 [10.63951, 13.14186]95
## 2 United States AAN     2017 N(12, 0.41)  12.0 [10.75667, 13.25656]95
resids <- augment(fit_us)
s <- sd(resids$.resid) 
mix_us %>%
  forecast(h=1)%>%
  mutate(low= .mean -1.96*s,
         high= .mean +1.96*s)
## # A fable: 2 x 7 [1Y]
## # Key:     Country, .model [2]
##   Country       .model  Year     Exports .mean   low  high
##   <fct>         <chr>  <dbl>      <dist> <dbl> <dbl> <dbl>
## 1 United States ANN     2017 N(12, 0.41)  11.9  10.7  13.1
## 2 United States AAN     2017 N(12, 0.41)  12.0  10.8  13.2

Again for both models the interval computed by R and manually are really close.

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.

Here is the plot of China’s GDP:

china_gdp <- global_economy %>%
  filter(Country=="China") %>%
  mutate(gdp_bill = GDP/1000000000)

china_gdp %>%
  autoplot(gdp_bill) + 
  labs(y = "GDP in billions", title = "China GDP in USD$"
       , x = "Date")

Let’s take transformation:

lamb_cna <- china_gdp %>%
  features(gdp_bill, features = guerrero) %>%
  pull(lambda_guerrero)

china_gdp %>%
  autoplot(box_cox(gdp_bill, lamb_cna)) + 
  labs(x = "Date"
       , y = "GDP Transformed from Billions of $"
       , title = "China GDP"
       , subtitle = "lambda = -0.03446284"
       ) +
  scale_y_continuous(labels = scales::number)

Let’s try our different ETS models and see how they are behaving

china_models <- china_gdp %>%
  model(
    ets1 = ETS(gdp_bill ~ error("A") + trend("A") + season("N")) # no seasonality
    , ets2 = ETS(gdp_bill ~ error("A") + trend("N") + season("N")) #remove the trend--mmm
    , ets_damp = ETS(gdp_bill ~ error("A")+trend("Ad")+season("N"))
    , ets_bc = ETS(box_cox(gdp_bill, lamb_cna))
    , ets_bc_damp = ETS(box_cox(gdp_bill, lamb_cna) ~ trend("Ad"))
    , ets_log = ETS(log(gdp_bill))
  )

The models are ready and let;s use them to forecast the next 10 years and plot the results:

china_models %>%
  forecast(h = 10) %>%
  autoplot(china_gdp, level = NULL) + 
  labs(x = "Date"
       , y = "GDP in Billions of $"
       , title = "China GDP"
       #, subtitle = "lambda = -0.03446284"
       ) +
  scale_y_continuous(labels = scales::dollar)

china_models %>% accuracy()
## # A tibble: 6 × 11
##   Country .model      .type      ME  RMSE   MAE   MPE  MAPE  MASE RMSSE     ACF1
##   <fct>   <chr>       <chr>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>
## 1 China   ets1        Traini…  23.6  190.  96.0 1.36   7.72 0.443 0.453  0.00882
## 2 China   ets2        Traini… 210.   416. 213.  8.25  10.9  0.983 0.991  0.789  
## 3 China   ets_damp    Traini…  29.4  190.  94.8 1.86   7.56 0.437 0.454 -0.00684
## 4 China   ets_bc      Traini… -44.5  301. 130.  0.451  7.24 0.601 0.719  0.677  
## 5 China   ets_bc_damp Traini… -18.5  221. 109.  1.17   7.23 0.503 0.528  0.346  
## 6 China   ets_log     Traini… -35.2  287. 125.  0.733  7.21 0.578 0.684  0.654

From the models above, when looking at RMSE, it appears that ets1, which is an AAN model is the best fit, along with ets_damp. This makes sense because it follows the trend. The box-cox transformation did not do much to help, although mixed with the trend (damp) model, although not great, do not appear as out of control as the ets, and log transformed model, which become exponential.

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?

Let’s filter and plot the data first:

gas_aus <- aus_production %>%
  filter(Quarter > yearquarter("1980 Q1")) %>%
  select(Quarter, Gas)
  
gas_aus %>% 
  autoplot(Gas) +
  labs(y = "Petajouls"
        , x = "Date"
        , title = "Gas Production in Australia")

let’s create models and plot the forecast for next 10 years:

models_gas <- 
  gas_aus %>%
  model(
    gas_m = ETS(Gas ~ error("M") + trend("A") + season("M"))
    , gas_damp = ETS(Gas ~ error("A") + trend("Ad") + season("M"))
    , gas_damp_sm = ETS(Gas ~ error("A") + trend("Ad",phi=0.90) + season("M"))
    )
 
gas_fc <-
  models_gas %>%
  forecast(h = 10)

plotting the forcast for all the models:

gas_fc %>% 
  autoplot(gas_aus, level = NULL) +
  labs(y = "Petajouls", title = "Gas Production in Australia"
       , x = "Date") +
  facet_grid(.model~.)

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

set.seed(54321)

aus_turnover <- aus_retail %>%
  filter(`Series ID` == sample(aus_retail$`Series ID`,1))

autoplot(aus_turnover, Turnover) +
  labs(title = "Turnover by Month")

a. Why is multiplicative seasonality necessary for this series?

Multiplicative seasonality is needed because the seasonal variations become larger as the trend changes. This method takes into account these changes–the variance–over time.

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

models_aus_turnover <- aus_turnover %>%
  model(
    multiplicative = ETS(Turnover ~ error("A") + trend("A") + season("M")),
    add_damped = ETS(Turnover ~ error("A") + trend("Ad") + season("M")),
    add_damped_0.95 = ETS(Turnover ~ error("A") + trend("Ad",phi=0.95) + season("M"))
    )
 
fc_aus_turnover <- models_aus_turnover %>%
  forecast(h = 60) 
  
(fc_aus_turnover %>%
  autoplot(aus_turnover, level = NULL) +
  labs(y = "Turnover", title = "Turnover by Month") +  facet_grid(.model~.))

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

test_models_aus_turnover<-
models_aus_turnover %>%
  accuracy() %>%
  select(.model,RMSE)

test_models_aus_turnover
## # A tibble: 3 × 2
##   .model           RMSE
##   <chr>           <dbl>
## 1 multiplicative   3.02
## 2 add_damped       2.98
## 3 add_damped_0.95  2.96

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

aus_turnover %>%
  model(damped = ETS(Turnover ~ error("A") + trend("Ad") + season("M"))) %>%
          gg_tsresiduals() + 
   labs(title = "Residuals, Damped Model")

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?

aus_turnover_train <- aus_turnover %>%
  filter(year(Month) < 2011)

model_train_aus_turnover <- aus_turnover_train %>%  # fit the new train set
  model(add_damped_train = ETS(Turnover ~ error("A") + trend("Ad") + season("M")))

## test the model

aus_turnover_test <- aus_turnover %>%
  filter(year(Month) > 2010)  #split to assess prior to these outliers

model_test_aus_turnover <- aus_turnover_test %>%  # fit the new train set
  model(add_damped_test = ETS(Turnover ~ error("A") + trend("Ad") + season("M")))

#VS excercise 7 from chapter 5.11

model_5.11 <- aus_turnover_train %>%
  model(SNAIVE(Turnover))

## Compare

aus_train_turnover <- model_train_aus_turnover %>%
  accuracy() %>%
  select(.model,RMSE)

aus_test_turnover <- model_test_aus_turnover %>%
  accuracy() %>%
  select(.model,RMSE)

mod_5.11 <- model_5.11 %>%
  accuracy() %>%
  select(.model,RMSE)
all <- rbind(aus_train_turnover,aus_test_turnover, mod_5.11)
all
## # A tibble: 3 × 2
##   .model            RMSE
##   <chr>            <dbl>
## 1 add_damped_train  2.46
## 2 add_damped_test   3.61
## 3 SNAIVE(Turnover)  5.64

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?

lambda_turnover <- 
  aus_turnover %>%
  features(Turnover, features = guerrero) %>%
  pull(lambda_guerrero)

aus_turnover_trans <- aus_turnover_train %>%
  mutate(bc = box_cox(Turnover,lambda_turnover))

stl_ets_holt <- aus_turnover_trans %>% 
  model("STL(BoxCox)" = STL(bc ~ season(window="periodic"),robust=T)
        ,"ETS (BoxCox)" = ETS(bc)
        )

#different response variable- separate

holts <- aus_turnover_trans %>%
    model("Holt Winters Damp" =ETS(Turnover ~ error("M") + trend("Ad") + season("M")))
      
all_n <- rbind(accuracy(stl_ets_holt),accuracy(holts))

all_n
## # A tibble: 3 × 12
##   State   Industry .model .type       ME   RMSE    MAE     MPE  MAPE  MASE RMSSE
##   <chr>   <chr>    <chr>  <chr>    <dbl>  <dbl>  <dbl>   <dbl> <dbl> <dbl> <dbl>
## 1 South … Clothin… STL(B… Trai… -1.78e-3 0.0427 0.0314 -0.0822  1.10 0.379 0.409
## 2 South … Clothin… ETS (… Trai… -6.37e-5 0.0511 0.0380 -0.0294  1.34 0.458 0.490
## 3 South … Clothin… Holt … Trai…  2.26e-1 2.51   1.70    0.113   5.30 0.428 0.445
## # ℹ 1 more variable: ACF1 <dbl>

Using the same data, and looking at RMSE, it appears that the transformed data performed considerably better