DATA624_HW5
William Aiken
10/6/2024
Exercise 8.1
Consider the the number of pigs slaughtered in Victoria, available in the aus_livestock dataset.
- 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.
The optimal values of alpha is 0.32 and the optimal value for lambda-zero is 100,646
pig_df <- aus_livestock %>% filter(State == 'Victoria' & Animal == 'Pigs')fit <- pig_df %>% model(ETS(Count ~ error("A") + trend("N") + season("N")))tidy(fit)## # A tibble: 2 × 5
## Animal State .model term estimate
## <fct> <fct> <chr> <chr> <dbl>
## 1 Pigs Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + sea… alpha 3.22e-1
## 2 Pigs Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + sea… l[0] 1.01e+5fc <- fit |>
forecast(h = 5)
kableExtra::kable(fc)| Animal | State | .model | Month | Count | .mean |
|---|---|---|---|---|---|
| Pigs | Victoria | ETS(Count ~ error(“A”) + trend(“N”) + season(“N”)) | 2019 Jan | N(95187, 8.7e+07) | 95186.56 |
| Pigs | Victoria | ETS(Count ~ error(“A”) + trend(“N”) + season(“N”)) | 2019 Feb | N(95187, 9.7e+07) | 95186.56 |
| Pigs | Victoria | ETS(Count ~ error(“A”) + trend(“N”) + season(“N”)) | 2019 Mar | N(95187, 1.1e+08) | 95186.56 |
| Pigs | Victoria | ETS(Count ~ error(“A”) + trend(“N”) + season(“N”)) | 2019 Apr | N(95187, 1.1e+08) | 95186.56 |
| Pigs | Victoria | ETS(Count ~ error(“A”) + trend(“N”) + season(“N”)) | 2019 May | N(95187, 1.2e+08) | 95186.56 |
- 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.
We can see that the interval produced by calculation is the same (some rounding error) generated by the forecast function.
temp <- fc$Count[1] %>% unlist()
temp[1] - temp[2]*1.96## mu
## 76854.45temp[1] + temp[2]*1.96## mu
## 113518.7hilo(fc)$`95%`[1]## <hilo[1]>
## [1] [76854.79, 113518.3]95Exercise 8.5
Data set global_economy contains the annual Exports from many countries. Select one country to analyse.
- Plot the Exports series and discuss the main features of the data.
There is an upward trend to the data with a sharp peak in 2000 and then a downward trend after that. There is no seasonality but potentially a cycle to the data.
canada_df <- global_economy %>% filter(Country == 'Canada')
canada_df %>% autoplot(Exports)
- Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
fit <- canada_df %>% model(ETS(Exports ~ error("A") + trend("N") + season("N")))
fit %>% forecast(h = 10) %>%
autoplot(canada_df)
- Compute the RMSE values for the training data.
The RMSE of the training data is 1.6
fit %>% 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 Canada "ETS(Exports ~ … Trai… 0.240 1.62 1.19 0.885 4.15 0.983 0.991 0.310- 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.
The RMSE for the ETS(A,A,N) model is slightly higher than the ETS(A,N,N) model. The first model is simpler but there appears to be a trend component to the data so I feel the second ETS(A,A,N) model may be more appropriate.
fit2 <- canada_df %>% model(ETS(Exports ~ error("A") + trend("A") + season("N")))
fit2 %>% 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 Canada "ETS(Exports… Trai… -0.0135 1.63 1.19 0.0538 4.22 0.981 0.996 0.160- Compare the forecasts from both methods. Which do you think is best?
Looking at the prediction of the second model, it takes the downward trend after 2000 into account. The first model is a naive model and assumes that the last value in the dataset is the best prediction for all future timepoints. I think that the second model may be more appropriate.
fit2 %>% forecast(h = 10) %>%
autoplot(canada_df)
- 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.
I’m getting tighter confidence intervals when I calculate them from the RMSE than the CI generated from the accuracy function.
fc1 <- canada_df %>% stretch_tsibble(.init = 57) %>% model(ETS(Exports ~ error("A") + trend("N") + season("N"))) %>% forecast(h=1)
acc1 <- canada_df %>% stretch_tsibble(.init = 57) %>% model(ETS(Exports ~ error("A") + trend("N") + season("N"))) %>% forecast(h=1) %>% accuracy(canada_df)## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing.
## 1 observation is missing at 2018temp1 <- acc1["RMSE"]
temp2 <- fc1[2, ".mean"]
names(temp2) <- "lower"
temp2 - temp1*1.96*sqrt(58)## lower
## 1 29.79085names(temp2) <- "upper"
temp2 + temp1*1.96*sqrt(58)## upper
## 1 31.99722hilo(fc1)$`95%`[2]## <hilo[1]>
## [1] [27.66725, 34.12082]95Now for the AAN model
fc2 <- canada_df %>% stretch_tsibble(.init = 57) %>% model(ETS(Exports ~ error("A") + trend("A") + season("N"))) %>% forecast(h=1)
acc2 <- canada_df %>% stretch_tsibble(.init = 57) %>% model(ETS(Exports ~ error("A") + trend("A") + season("N"))) %>% forecast(h=1) %>% accuracy(canada_df)## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing.
## 1 observation is missing at 2018temp1 <- acc2["RMSE"]
temp2 <- fc2[2, ".mean"]
names(temp2) <- "lower"
temp2 - temp1*1.96*sqrt(58)## lower
## 1 29.98082names(temp2) <- "upper"
temp2 + temp1*1.96*sqrt(58)## upper
## 1 31.57869hilo(fc2)$`95%`[2]## <hilo[1]>
## [1] [27.47781, 34.0817]95Exercise 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.]
To start we can see that China has experienced exponential growth in their GDP since the 1980s with it really taking off in the 90s.
chinese_df <- global_economy %>% filter(Country == 'China') %>% select(Country, Year, GDP)
chinese_df %>% autoplot()## Plot variable not specified, automatically selected `.vars = GDP`
Using the AAN (Holt) model we get a very optimistic forecast for China’s future GDP doubling in the next 50 years.
fit <- chinese_df %>% model(ETS(GDP ~ error("A") + trend("A") + season("N"))) %>% forecast(h=50)
fit %>% autoplot(chinese_df)
The damped model gives us a more realistic model with growth that tappers off with time.
fit <- chinese_df %>% model(ETS(GDP ~ error("A") + trend("Ad") + season("N"))) %>% forecast(h=50)
fit %>% autoplot(chinese_df)
The Box-Cox transformation makes the Chinese GDP plot linear, with the AAN (Holt) model continuing this linear trend
library(caret)## Loading required package: lattice##
## Attaching package: 'caret'## The following objects are masked from 'package:fabletools':
##
## MAE, RMSEchinese_trans <- chinese_df %>% select(Year, GDP) %>% as.data.frame()
bc_trans <- preProcess(chinese_trans["GDP"], method = c("BoxCox"))
transformed <- predict(bc_trans, chinese_trans["GDP"])
chinese_trans$GDP_t <- transformed$GDP
chinese_trans <- as_tsibble(chinese_trans, index = 'Year')
fit <- chinese_trans %>% model(ETS(GDP_t ~ error("A") + trend("A") + season("N"))) %>% forecast(h=50)
fit %>% autoplot(chinese_trans)
Lastly, let’s plot the Box-Cox transformed Chinese GDP data with a damped model. The Box-Cox transformation does really change the shape of the prediction, just the scale.
fit <- chinese_trans %>% model(ETS(GDP_t ~ error("A") + trend("Ad") + season("N"))) %>% forecast(h=50)
fit %>% autoplot(chinese_trans)
### Exercise 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?
By looking at the plot of the time series and prior assignments, we see that the the seasonality is not constant across the series. There is heteroskedasticity in the variance of the seasonality, this is going to require a multiplicative seasonality.
gas_hw <- aus_production %>% select(Quarter, Gas)
gas_hw %>% autoplot()## Plot variable not specified, automatically selected `.vars = Gas`
We are going to use the ‘model’ function to select an optimal model for us. We see that the optimal model has multiplicative error, additive error and multiplicative seasonality. I’m using a large h in my forecasts so that I can compare these forecasts to the dampened forecasts later on.
fit <- gas_hw %>% model(ETS(Gas))
report(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.389fc <- gas_hw %>% model(ETS(Gas ~ error("M") + trend("A") + season("M"))) %>% forecast(h=50)
fc %>% autoplot(gas_hw)
Lastly, let’s replace the ‘additive’ trend with a dampened additive trend. The trend of the dampened forecasts are slightly more conservative than the optimal model.
fc <- gas_hw %>% model(ETS(Gas ~ error("M") + trend("Ad") + season("M"))) %>% forecast(h=50)
fc %>% autoplot(gas_hw)
Exercise 8.8
Recall your retail time series data (from Exercise 7 in Section 2.10).
- Why is multiplicative seasonality necessary for this series?
Similar to the aus_production Gas data set, there is uneven variance is the seasonality of the time series requiring multiplicative seasonality
set.seed(555)
myseries <- aus_retail |>
filter(`Series ID` == sample(aus_retail$`Series ID`,1))
myseries |> autoplot()## Plot variable not specified, automatically selected `.vars = Turnover`
- Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
fc <- myseries %>% model(ETS(Turnover ~ error("M") + trend("A") + season("M"))) %>% forecast(h=50)
fc %>% autoplot(myseries)
Let’s dampened the forecast, very similar trend
fc_d <- myseries %>% model(ETS(Turnover ~ error("M") + trend("Ad") + season("M"))) %>% forecast(h=50)
fc_d %>% autoplot(myseries)
- Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer? The RMSE of the dampened Multiplicative Holt-Winter is smaller so I would select that method.
fit <- myseries %>% model( HW_damp = ETS(Turnover ~ error("M") + trend("Ad") + season("M")), HW = ETS(Turnover ~ error("M") + trend("A") + season("M")))
kableExtra::kable(accuracy(fit) %>% select(c('.model','RMSE')))| .model | RMSE |
|---|---|
| HW_damp | 5.529967 |
| HW | 5.532134 |
- Check that the residuals from the best method look like white noise. The residuals look evenly distributed, and homoscedastic. There is some autocorrelation which I believe is to be expected with seasonal data.
fit %>% select(HW) %>% gg_tsresiduals()
- 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?
train <- myseries %>% filter(year(Month) < '2011') %>%
model(naive = SNAIVE(Turnover),
HW = ETS(Turnover ~ error("M") + trend("A") + season("M"))
)
fc_train <- train %>%
forecast(h=120)
fc_train %>% autoplot(myseries)
These results are a little surprising, the RMSE of the SNAIVE model is much small than that of the Multiplicative Holt Winter. I think that a dampened muliplicative Holt-Winter method would have performed better.
kableExtra::kable(accuracy(fc_train, myseries) %>% select(c('.model','RMSE')))## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing.
## 24 observations are missing between 2019 Jan and 2020 Dec| .model | RMSE |
|---|---|
| HW | 23.203162 |
| naive | 7.965871 |
Exercise 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?
- Note to self, the ‘preProcess’ function doesn’t like tsibble objects. They need to be cast to df first
The level is linear after the Box-Cox transformation. The slope is a tiny component and the seasonal component becomes less peaked over time.
library(caret)
myseries_df <- as.data.frame(myseries)
bc_trans <- preProcess(myseries_df["Turnover"], method = c("BoxCox"))
transformed <- predict(bc_trans, myseries_df["Turnover"])
myseries_df$Turnover_t <- transformed$Turnover
myseries_df <- as_tsibble(myseries_df, index = 'Month')
myseries_df %>% filter(year(Month) < '2011') %>% model(ETS(Turnover_t ~ error("M") + trend("A") + season("M"))) %>%
components() |>
autoplot() +
labs(title = "Classical additive decomposition of retail turnover")## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).
fit <- myseries_df %>% filter(year(Month) < '2011') %>% model(BC = ETS(Turnover_t ~ error("M") + trend("A") + season("M"))) %>% forecast(h=120)
fit %>% autoplot(myseries_df)
The RMSE is shrunk considerably. I believe that the RMSE of the transformed model can be compared directly to the RMSE of the prior models
kableExtra::kable(accuracy(fit, myseries_df) %>% select(c('.model','RMSE')))## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing.
## 24 observations are missing between 2019 Jan and 2020 Dec| .model | RMSE |
|---|---|
| BC | 0.2614963 |