624 Homework Five
Raymond Fries
2023-02-28
library(fpp3)## ── Attaching packages ────────────────────────────────────────────── fpp3 0.5 ──## ✔ tibble 3.1.8 ✔ tsibble 1.1.3
## ✔ dplyr 1.1.0 ✔ tsibbledata 0.4.1
## ✔ tidyr 1.3.0 ✔ feasts 0.3.0
## ✔ lubridate 1.9.2 ✔ fable 0.3.2
## ✔ ggplot2 3.4.1 ✔ fabletools 0.3.2## ── 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()Exercise 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.
The ETS model will choose the best model, so if you don’t include trend or season in the model, it will still try those components and give back the most optimal model. By default the model uses log-likelihood as the optimization metric, although it can be set to mse, amse, sigma or mae.
Using the Australian livestock dataset the ETS model gave an optimal setting of error A(additive), trend N(None) and seasonal A(additive). The optimal alpha was .3579401, while the optimal ℓ0 was 95487.5.
data(aus_livestock)
unique(aus_livestock$Animal)## [1] Bulls, bullocks and steers Calves
## [3] Cattle (excl. calves) Cows and heifers
## [5] Lambs Pigs
## [7] Sheep
## 7 Levels: Bulls, bullocks and steers Calves ... Sheepvic_pigs <- aus_livestock %>% filter(State=='Victoria',Animal=='Pigs')fit <- vic_pigs |>
model(ETS(Count ~ error() ))
fc <- fit |>
forecast(h = 4)
fc %>% autoplot(vic_pigs) +
geom_line(aes(y = .fitted), col="orange",
data = augment(fit)) +
labs(y="Number of Pigs Slaughtered", title="Victoria Slaughtered Pigs") +
guides(colour = "none")
autoplot(fc)
report(fit)## Series: Count
## Model: ETS(A,N,A)
## Smoothing parameters:
## alpha = 0.3579401
## gamma = 0.0001000139
##
## Initial states:
## l[0] s[0] s[-1] s[-2] s[-3] s[-4] s[-5] s[-6]
## 95487.5 2066.235 6717.311 -2809.562 -1341.299 -7750.173 -8483.418 5610.89
## s[-7] s[-8] s[-9] s[-10] s[-11]
## -579.8107 1215.464 -2827.091 1739.465 6441.989
##
## sigma^2: 60742898
##
## AIC AICc BIC
## 13545.38 13546.26 13610.24b) Compute a 95% prediction interval for the first forecast using y ± 1.96 s, where s is the standard deviation of the residuals. Compare your interval with the interval produced by R.
The difference between the interval I did and the interval the forecast did is +- .727
std <- 7794
mu <- fc$.mean[[1]][1]
lower <- mu - 1.96 * std
upper <- mu + 1.96 * stdpaste0("Lower Bound: ",round(lower,3),' Upper Bound: ',round(upper,3))## [1] "Lower Bound: 69148.467 Upper Bound: 99700.947"nf <- unpack_hilo(hilo(fc, 95))[1,]
nf[,c("Month","Count",".mean",'95%')]## # A tsibble: 1 x 4 [1M]
## Month Count .mean `95%`
## <mth> <dist> <dbl> <hilo>
## 1 2019 Jan N(84425, 6.1e+07) 84425. [69149.19, 99700.22]95Exercise 8.5
Data set global_economy contains the annual Exports from many countries. Select one country to analyse.
data(global_economy)a) Plot the Exports series and discuss the main features of the data.
The main feature of the data is an upward trend. There does not appear to be any seasonality since the peaks and valleys all appear to have different lengths in time. Also, the US Exports is missing the last period, which has to be removed in order to forecast.
us_exports <- global_economy %>%
filter(Country == "United States")
summary(us_exports)## Country Code Year GDP
## United States :58 USA :58 Min. :1960 Min. :5.433e+11
## Afghanistan : 0 ABW : 0 1st Qu.:1974 1st Qu.:1.584e+12
## Albania : 0 AFG : 0 Median :1988 Median :5.455e+12
## Algeria : 0 AGO : 0 Mean :1988 Mean :6.992e+12
## American Samoa: 0 ALB : 0 3rd Qu.:2003 3rd Qu.:1.138e+13
## Andorra : 0 AND : 0 Max. :2017 Max. :1.939e+13
## (Other) : 0 (Other): 0
## Growth CPI Imports Exports
## Min. :-2.776 Min. : 13.56 Min. : 4.030 Min. : 4.809
## 1st Qu.: 1.919 1st Qu.: 23.13 1st Qu.: 7.267 1st Qu.: 6.976
## Median : 3.176 Median : 55.54 Median :10.265 Median : 9.038
## Mean : 3.043 Mean : 56.85 Mean :10.378 Mean : 8.763
## 3rd Qu.: 4.450 3rd Qu.: 83.89 3rd Qu.:13.413 3rd Qu.:10.606
## Max. : 7.259 Max. :112.41 Max. :17.427 Max. :13.639
## NA's :1 NA's :1 NA's :1
## Population
## Min. :180671000
## 1st Qu.:214383750
## Median :245659000
## Mean :251274672
## 3rd Qu.:289487248
## Max. :325719178
## us_exports$Exports## [1] 4.969630 4.899698 4.809122 4.870028 5.103529 4.988571 5.018405
## [8] 5.048161 5.082228 5.088734 5.549855 5.391815 5.524118 6.669005
## [15] 8.177231 8.212749 7.963146 7.639201 7.930379 8.743028 9.808647
## [22] 9.506172 8.466660 7.613677 7.483444 6.975559 6.993228 7.472850
## [29] 8.464352 8.913368 9.229297 9.636020 9.680747 9.519216 9.864047
## [36] 10.605515 10.710722 11.079797 10.484799 10.268281 10.664643 9.666071
## [43] 9.132386 9.037519 9.625368 9.996398 10.654792 11.497907 12.514398
## [50] 11.011656 12.378301 13.573792 13.606608 13.639312 13.620044 12.499044
## [57] 11.890622 NAus_exports %>% autoplot(Exports)
b) Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
us_fit <- head(us_exports,-1) |>
model(ETS(Exports ~ error("A") + trend("N") + season("N")))
us_fc <- us_fit |>
forecast(h = 4)us_fc %>% autoplot(us_exports) +
geom_line(aes(y = .fitted), col="orange",
data = augment(us_fit)) +
labs(title="US Exports") +
guides(colour = "none")
c) Compute the RMSE values for the training data.
accuracy(us_fit)$RMSE## [1] 0.6270672d) 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 results from the two RMSEs are basically the same, only differing by .0121, in favor of the AAN ets model. The biggest difference is that there is clearly an upward trend in the data set. With the forecast for the ANN ets model the mean of the forecast is flat, while the mean for the AAN model is trends upwards. It’s basically the same as naive versus rw drift.
us_fit2 <- head(us_exports,-1) |>
model(ETS(Exports ~ error("A") + trend("A") + season("N")))
us_fc2 <- us_fit2 |>
forecast(h = 4)us_fc2 %>% autoplot(us_exports) +
geom_line(aes(y = .fitted), col="orange",
data = augment(us_fit2)) +
labs(title="US Exports") +
guides(colour = "none")
accuracy(us_fit2)$RMSE## [1] 0.6149566e) Compare the forecasts from both methods. Which do you think is best?
Long term, I think the AAN model is best as it accounts for the upward trend in the data set. It would be unreasonable to think that the US Exports mean is going to flatline.
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.
The intervals that I calculated were off by roughly .036 for the first model and roughly .035 for the second model.
mu1 <- us_fc$.mean[1]
sd1 <- accuracy(us_fit)$RMSE
lower_bound <- mu1 - sd1*1.96
upper_bound <- mu1 + sd1*1.96
mu2 <- us_fc2$.mean[1]
sd2 <- accuracy(us_fit2)$RMSE
lower_bound2 <- mu2 - sd2*1.96
upper_bound2 <- mu2 + sd2*1.96paste0("Model 1: Lower Bound: ",round(lower_bound,3),' Upper Bound: ',round(upper_bound,3))## [1] "Model 1: Lower Bound: 10.662 Upper Bound: 13.12"paste0("Model 1: Lower Bound: ",round(lower_bound2,3),' Upper Bound: ',round(upper_bound2,3))## [1] "Model 1: Lower Bound: 10.801 Upper Bound: 13.212"us_fc %>% hilo(95) %>% pull('95%') %>% head(1)## <hilo[1]>
## [1] [10.63951, 13.14186]95us_fc2 %>% hilo(95) %>% pull('95%') %>% head(1)## <hilo[1]>
## [1] [10.75667, 13.25656]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.]
china <- global_economy %>% filter(Country=="China")
china %>% autoplot(GDP)
lambda <- global_economy %>%
filter(Country == "China") %>%
features(GDP, features = guerrero) %>%
pull(lambda_guerrero)china_fit <- china %>%
model(`normal` = ETS(GDP ~ error('A') + trend('N') + season('N')),
`bc` = ETS(box_cox(GDP,lambda) ~ error("A") + trend("A") + season("N")),
`bc_damped` = ETS(box_cox(GDP,lambda) ~ error("A") + trend("Ad", phi = 0.90) + season("N")),
`holts_method` = ETS(GDP ~ error('A') + trend('A') + season('N')),
`damped_Holts_method` = ETS(GDP ~ error('A') + trend('Ad', phi = 0.9) + season('N')),
)
china_fc <- china_fit %>%
forecast(h = 20)
china_fc %>%
autoplot(global_economy, level = NULL) +
labs(title="China GDP Forecasts") 
Since there wasn’t a lot of variance in the China GDP data, the box-cox was able to make a pretty nice parabola. The damping effect on the box-cox with a phi of .9, has a nice projection with a slight upward curve. The normal ets is flat, or possibly even declining. Holt’s methods both have decent projections, with the dampening effect clearly restricting the normal Holt method.
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?
Multiplicative seasonality is necessary because there is the pattern of what looks like two quarters up to two quarters down and the variance is growing as the time series goes on. Additive would not capture that growth in variance, so the additive would underestimate future quarters.
Dampening the data just had the effect of leveling out the arc of the highs and lows.So it basically removed growth in the variance as the forecast went further out, as the variance of the damped forecast looks the same as in the regular multiplicative model.
aus_gas <- aus_production %>% select(Gas)
gas_fit <- aus_gas %>%
model(initial_model = ETS(Gas ~ error() + trend() + season()))
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.389gas_fit %>%
forecast(h=8) %>%
autoplot(aus_production, level = NULL) +
labs(title="Australian Gas Production")
damped_gas <- aus_production %>%
model(damped_mult = ETS(Gas ~ error("M") + trend("Ad", phi = 0.9) + season("M")))
damped_gas_fit <- aus_gas %>% model(initial_model = ETS(Gas ~ error('M') + trend('Ad',phi=.9) + season('M')))
report(damped_gas_fit)## Series: Gas
## Model: ETS(M,Ad,M)
## Smoothing parameters:
## alpha = 0.5672292
## beta = 0.2225463
## gamma = 0.09061058
## phi = 0.9
##
## Initial states:
## l[0] b[0] s[0] s[-1] s[-2] s[-3]
## 5.75995 0.09449598 0.929894 1.174977 1.074919 0.82021
##
## sigma^2: 0.0034
##
## AIC AICc BIC
## 1688.056 1688.921 1718.516damped_gas_fit %>%
forecast(h=8) %>%
autoplot(aus_production, level = NULL) +
labs(title="Australian Gas Production")
Exercise 8.8
Recall your retail time series data (from Exercise 8 in Section 2.10).
a) Why is multiplicative seasonality necessary for this series?
Multiplicative seasonality is necessary in this series because the variance isn’t stationary and grows and shrinks at different rates throughout the timeseries.
set.seed(12345678)
myseries <- aus_retail |>
filter(`Series ID` == sample(aus_retail$`Series ID`,1))
myseries %>% autoplot(Turnover)
myseries_train <- myseries %>%
filter(year(Month) < 2017)
myseries_train %>% autoplot(Turnover)
b) Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
The dampening effect with a phi value of .2 compared to a phi setting of point .9 is only different in the aspect that the lower forecast with phi of .2 has a slightly lower bottom than .9. Both peaks appear to be the same.
myseries_fit <- myseries_train %>%
model(ETS(Turnover ~ error() + trend() + season()))
myseries_fit %>%
forecast(h = 13) %>%
autoplot(myseries,level=NULL,color='blue')+
labs(title="Holt Winters auto chosen phi=0.9794538 ")
hw_fit <- myseries_train %>% model(ETS(Turnover ~ error('M') + trend('Ad') + season('M')))
report(hw_fit)## Series: Turnover
## Model: ETS(M,Ad,M)
## Smoothing parameters:
## alpha = 0.509344
## beta = 0.0001129272
## gamma = 0.100066
## phi = 0.9797018
##
## Initial states:
## l[0] b[0] s[0] s[-1] s[-2] s[-3] s[-4] s[-5]
## 2.370523 0.1042923 0.8065457 0.7613085 0.8265112 1.337986 0.9988949 1.0556
## s[-6] s[-7] s[-8] s[-9] s[-10] s[-11]
## 1.088849 1.103015 1.13785 1.0171 1.000835 0.8655057
##
## sigma^2: 0.0048
##
## AIC AICc BIC
## 1636.973 1639.071 1706.157hw_fit %>%
forecast(h = 13) %>%
autoplot(myseries,level=NULL)+
labs(title="Holt Winters auto chosen phi=0.9794538 ")
hwd2_fit <- myseries_train %>% model(ETS(Turnover ~ error('M') + trend('Ad',phi=.2) + season('M')))
hwd2_fit %>%
forecast(h = 13) %>%
autoplot(myseries,level=NULL)+
labs(title="Holt Winters phi=.2")
hwd9_fit <- myseries_train %>% model(ETS(Turnover ~ error('M') + trend('Ad',phi=.6) + season('M')))
hwd9_fit %>%
forecast(h = 13) %>%
autoplot(myseries,level=NULL) +
labs(title="Holt Winters phi=.6")
c) Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
The Holt-Winters method has a slightly better rmse score at .5223 compared to .5263. Looking at the plots the Holt-Winters method with phi=.97 is definitely better as it matches the test section almost perfectly.
msf <-accuracy(myseries_fit)
hwf <- accuracy(hw_fit)paste0("Myseries RMSE ",round(msf[,11],4),' Holt Winters phi=.97 RMSE ',round(hwf[,11],4))## [1] "Myseries RMSE 0.5213 Holt Winters phi=.97 RMSE 0.5279"d) Check that the residuals from the best method look like white noise.
Yes, the residuals look like white noise and the distribution of the residuals looks normal.
hw_fit %>%
gg_tsresiduals() +
ggtitle("Multiplicative Method")
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?
The Holt-Winters method is definitely better than the seasonal naive method. As the rmse of the Holt-Winters method is 1.15 and the seasonal naive’s is 1.55. Also on the plot’s we can see that the Holt-Winters forecast traces the test data more closely.
train2 <- myseries %>%
filter(year(Month) < 2011)
fit_train <- train2 %>%
model(hw = ETS(Turnover ~ error('M') + trend('Ad') + season('M')),
snaive = SNAIVE(Turnover))
fc2 <- fit_train %>%
forecast(new_data = anti_join(myseries, train2))## Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover)`fc2 %>% autoplot(myseries, level = NULL)
fc2 %>% accuracy(myseries) %>%
select(.type, .model, RMSE)## # A tibble: 2 × 3
## .type .model RMSE
## <chr> <chr> <dbl>
## 1 Test hw 1.15
## 2 Test snaive 1.55Exercise 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?
The rmse is considerably lower for the stl decomposition with box-cox data, but that doesn’t really mean much because the rmse is lower because the values of the transformed data are lower. When comparing the plots of the Holt-Winters forecast and comparing it to the stl decomposition forecast plots, there is a clear difference in how the forecasts trace the test data. So the lower rmse does not mean a better model, what it does mean is that the data is on a different scale.
myseries2 <- myseries
lambda <- myseries2 %>%
features(Turnover, features = guerrero) %>%
pull(lambda_guerrero)
myseries2 %>%
model(STL(box_cox(Turnover,lambda) ~ season(window = "periodic"), robust = TRUE)) %>%
components() %>%
autoplot() 
decomp <- myseries2 %>%
model(STL_box = STL(box_cox(Turnover,lambda) ~ season(window = "periodic"), robust = TRUE)) %>%
components()
myseries2$Turnover <- decomp$season_adjusttrain2 <- myseries2 %>% filter(year(Month) < 2011)
fit3 <- train2 %>%
model(ETS(Turnover ~ error("M") + trend("Ad", phi=0.979996) + season("M")))fit3 %>% gg_tsresiduals() +
ggtitle("Residual Plot")
fc3 <- fit3 %>%
forecast(new_data = anti_join(myseries2, train2))## Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover)`fc3 %>% autoplot(myseries2, level = NULL)
fc3 %>% accuracy(myseries2) %>%
select(.model, .type, RMSE)## # A tibble: 1 × 3
## .model .type RMSE
## <chr> <chr> <dbl>
## 1 "ETS(Turnover ~ error(\"M\") + trend(\"Ad\", phi = 0.979996) + se… Test 0.117fc2 %>% accuracy(myseries) %>%
select(.type, .model, RMSE)## # A tibble: 2 × 3
## .type .model RMSE
## <chr> <chr> <dbl>
## 1 Test hw 1.15
## 2 Test snaive 1.55