DATA624_HW5
Gabriel Santos
2023-03-03
library(fpp2)
library(fpp3)
library(forecast)
library(ggplot2)
library(dplyr)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 ℓ, and generate forecasts for the next four months.
data("pigs")## Warning in data("pigs"): data set 'pigs' not foundhead((pigs))## Jan Feb Mar Apr May Jun
## 1980 76378 71947 33873 96428 105084 95741summary(pigs)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 33873 79080 91662 90640 101493 120184paste0("Frequency: ", frequency(pigs))## [1] "Frequency: 12"autoplot(pigs, colour = 'blue')
ggseasonplot(pigs)
ggsubseriesplot(pigs)
gglagplot(pigs)
Acf(pigs, lag.max=150)
Reviewing the graphs, he observed some cyclicality. I don’t see any
trend or seasonality in the data set. I’m going to use the SES (simple
exponential smoothing) function to explore the data set.
ses.pigs <- ses(pigs)
summary(ses.pigs)##
## Forecast method: Simple exponential smoothing
##
## Model Information:
## Simple exponential smoothing
##
## Call:
## ses(y = pigs)
##
## Smoothing parameters:
## alpha = 0.2971
##
## Initial states:
## l = 77260.0561
##
## sigma: 10308.58
##
## AIC AICc BIC
## 4462.955 4463.086 4472.665
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 385.8721 10253.6 7961.383 -0.922652 9.274016 0.7966249 0.01282239
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Sep 1995 98816.41 85605.43 112027.4 78611.97 119020.8
## Oct 1995 98816.41 85034.52 112598.3 77738.83 119894.0
## Nov 1995 98816.41 84486.34 113146.5 76900.46 120732.4
## Dec 1995 98816.41 83958.37 113674.4 76092.99 121539.8
## Jan 1996 98816.41 83448.52 114184.3 75313.25 122319.6
## Feb 1996 98816.41 82955.06 114677.8 74558.56 123074.2
## Mar 1996 98816.41 82476.49 115156.3 73826.66 123806.2
## Apr 1996 98816.41 82011.54 115621.3 73115.58 124517.2
## May 1996 98816.41 81559.12 116073.7 72423.66 125209.2
## Jun 1996 98816.41 81118.26 116514.6 71749.42 125883.4SES function: Smoothing parameters: alpha = 0.2971
ses_pigs.4 <- ses(pigs, h = 4)
ses_pigs.4$model## Simple exponential smoothing
##
## Call:
## ses(y = pigs, h = 4)
##
## Smoothing parameters:
## alpha = 0.2971
##
## Initial states:
## l = 77260.0561
##
## sigma: 10308.58
##
## AIC AICc BIC
## 4462.955 4463.086 4472.665- 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.
autoplot(ses.pigs) +
autolayer(fitted(ses.pigs), series="Fitted")
s <- sd(ses.pigs$residuals)
(ses.pigs$mean[1] + 1.96*s)## [1] 118952.8(ses.pigs$mean[1] - 1.96*s)## [1] 78679.97The results were slightly different than the results with the SES function.
ses.pigs## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Sep 1995 98816.41 85605.43 112027.4 78611.97 119020.8
## Oct 1995 98816.41 85034.52 112598.3 77738.83 119894.0
## Nov 1995 98816.41 84486.34 113146.5 76900.46 120732.4
## Dec 1995 98816.41 83958.37 113674.4 76092.99 121539.8
## Jan 1996 98816.41 83448.52 114184.3 75313.25 122319.6
## Feb 1996 98816.41 82955.06 114677.8 74558.56 123074.2
## Mar 1996 98816.41 82476.49 115156.3 73826.66 123806.2
## Apr 1996 98816.41 82011.54 115621.3 73115.58 124517.2
## May 1996 98816.41 81559.12 116073.7 72423.66 125209.2
## Jun 1996 98816.41 81118.26 116514.6 71749.42 125883.4s <- sd(ses.pigs$residuals)
ses.pigs$mean[1] + 1.96 * s ## [1] 118952.895% prediction interval computation
ses.pigs$mean[1] - 1.96 * s## [1] 78679.97Between the lower and upper limits there is a slightly small difference. Comparing my interval with the interval produced by R, they are similar.
Exercise 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.
fra_exports <- global_economy %>%
filter(Country == 'France')
head(fra_exports)## # 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 France FRA 1960 62651474947. NA 10.4 12.6 14.4 46814237
## 2 France FRA 1961 68346741504. 5.51 10.7 12.3 13.9 47444751
## 3 France FRA 1962 76313782252. 6.67 11.3 12.1 12.8 48119649
## 4 France FRA 1963 85551113767. 5.35 11.8 12.5 12.6 48803680
## 5 France FRA 1964 94906593388. 6.52 12.2 13.0 12.6 49449403
## 6 France FRA 1965 102160571409. 4.78 12.5 12.6 13.2 50023774Exports of goods and services (% of GDP)
fra_exports %>%
autoplot(Exports, colour ='blue') +
labs(y="% GDP", title="Exports France")
b.Use an ETS(A,N,N) model to forecast the series, and plot the forecasts
fit <- fra_exports %>%
model(ETS(Exports ~ error("A") + trend("N") + season("N")))
report(fit)## Series: Exports
## Model: ETS(A,N,N)
## Smoothing parameters:
## alpha = 0.9998995
##
## Initial states:
## l[0]
## 14.3989
##
## sigma^2: 1.3745
##
## AIC AICc BIC
## 257.9200 258.3644 264.1013- Compute the RMSE values for the training data.
fc <- fit %>%
forecast(h = 10)
fc %>%
autoplot(fra_exports, colour = 'blue') +
labs(y="% GDP", title="Exports France") +
guides(colour = "none")
head(fc)## # A fable: 6 x 5 [1Y]
## # Key: Country, .model [1]
## Country .model Year Exports .mean
## <fct> <chr> <dbl> <dist> <dbl>
## 1 France "ETS(Exports ~ error(\"A\") + trend(\"N\") + s… 2018 N(31, 1.4) 30.9
## 2 France "ETS(Exports ~ error(\"A\") + trend(\"N\") + s… 2019 N(31, 2.7) 30.9
## 3 France "ETS(Exports ~ error(\"A\") + trend(\"N\") + s… 2020 N(31, 4.1) 30.9
## 4 France "ETS(Exports ~ error(\"A\") + trend(\"N\") + s… 2021 N(31, 5.5) 30.9
## 5 France "ETS(Exports ~ error(\"A\") + trend(\"N\") + s… 2022 N(31, 6.9) 30.9
## 6 France "ETS(Exports ~ error(\"A\") + trend(\"N\") + s… 2023 N(31, 8.2) 30.9fit_acc <- accuracy(fit)
fit_acc## # 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 France "ETS(Exports… Trai… 0.284 1.15 0.840 1.17 3.86 0.983 0.991 -0.00542fit_2 <- fra_exports %>%
model(ETS(Exports ~ error("A") + trend("A") + season("N")))
fit_2_acc <- accuracy(fit_2)
fit_2_acc## # 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 France "ETS(Export… Trai… -0.0111 1.12 0.800 -0.243 3.81 0.936 0.964 0.0264d.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.
fc_2 <- fit_2 %>%
forecast(h = 10)
fc_2 %>%
autoplot(fra_exports, colour = 'blue') +
labs(y="% GDP", title="Exports France") +
guides(colour = "none")
e.Compare the forecasts from both methods. Which do you think is
best?
fc_2## # A fable: 10 x 5 [1Y]
## # Key: Country, .model [1]
## Country .model Year Exports .mean
## <fct> <chr> <dbl> <dist> <dbl>
## 1 France "ETS(Exports ~ error(\"A\") + trend(\"A\") + … 2018 N(31, 1.3) 31.2
## 2 France "ETS(Exports ~ error(\"A\") + trend(\"A\") + … 2019 N(31, 2.6) 31.5
## 3 France "ETS(Exports ~ error(\"A\") + trend(\"A\") + … 2020 N(32, 3.8) 31.8
## 4 France "ETS(Exports ~ error(\"A\") + trend(\"A\") + … 2021 N(32, 5) 32.1
## 5 France "ETS(Exports ~ error(\"A\") + trend(\"A\") + … 2022 N(32, 6.2) 32.4
## 6 France "ETS(Exports ~ error(\"A\") + trend(\"A\") + … 2023 N(33, 7.4) 32.7
## 7 France "ETS(Exports ~ error(\"A\") + trend(\"A\") + … 2024 N(33, 8.6) 33.0
## 8 France "ETS(Exports ~ error(\"A\") + trend(\"A\") + … 2025 N(33, 9.8) 33.4
## 9 France "ETS(Exports ~ error(\"A\") + trend(\"A\") + … 2026 N(34, 11) 33.7
## 10 France "ETS(Exports ~ error(\"A\") + trend(\"A\") + … 2027 N(34, 12) 34.0The predicted values for ETS(A,A,N) increase due to the additive trend parameter. It has a constant increase and is in line with the growing trend of export data. I think the ETS (A,A,N) model is adequate for forecasting this data set.
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.
95% prediction interval for the first forecast: ETS(A,N,N)
y_hat <- fc$.mean[1]
lower_limit_95_fc <- y_hat - (fit_acc$RMSE * 1.96)
upper_limit_95_fc <- y_hat + (fit_acc$RMSE * 1.96)
ets_ann_interval <- c(lower_limit_95_fc, upper_limit_95_fc)
ets_ann_interval## [1] 28.62385 33.1397195% prediction interval for the second forecast: ETS(A,A,N)
y_hat_2 <- fc_2$.mean[1]
lower_limit_95_fc2 <- y_hat_2 - (fit_2_acc$RMSE * 1.96)
upper_limit_95_fc2 <- y_hat_2 + (fit_2_acc$RMSE * 1.96)
ets_aan_interval <- c(lower_limit_95_fc2, upper_limit_95_fc2)
ets_aan_interval## [1] 28.97896 33.36920The model forecast 95% interval
fc_hilo <- fc %>% hilo()
fc_hilo$`95%`[1]## <hilo[1]>
## [1] [28.58393, 33.17963]95The R produced values are 0.04 lower for the lower limit and 0.04 higher for the higher limit.
The model forecast 95% intervals
fc_2_hilo <- fc_2 %>% hilo()
fc_2_hilo$`95%`[1]## <hilo[1]>
## [1] [28.89915, 33.44901]95The R produced values are 0.08 lower for the lower limit and 0.08 higher for the higher limit.
The confidence intervals produced by R seem a bit wider than the calculated values. The differences are very small, the confidence intervals calculated and produced by R seem to be similar to each other.
Exercise 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.]
options(scipen = 999)
chineseGDP <- global_economy %>%
filter(Country == 'China')
chineseGDP %>% autoplot(GDP, colour = 'blue') +
labs(title = 'Chinese GDP')
Get the optimized lambda value for BoxCox transformations.
lambda <- chineseGDP %>%
features(GDP, features = guerrero) %>%
pull(lambda_guerrero)
chineseGDPEtsOptionsComparision <- chineseGDP %>%
model(
ETS = ETS(GDP),
ETSBoxCox = ETS(box_cox(GDP, lambda)),
ETSDamped = ETS(GDP ~ trend('Ad', phi = 0.9)),
ETSLog = ETS(log(GDP))
)
chineseGDPEtsOptionsComparision %>%
forecast(h = 20) %>%
autoplot(chineseGDP, level = NULL) +
labs(title = 'ETS Forecast options comparison', subtitle = 'Chinese GDP')
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?
aus_production %>%
autoplot(Gas, colour = 'blue')
ETS model
fit <- aus_production %>%
model(fit = 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.389fit %>%
forecast(h = 4) %>%
autoplot(aus_production, colour = 'blue')
Exercise 8.8
Recall your retail time series data (from Exercise 8 in Section 2.10).
set.seed(1975)
myseries <- aus_retail %>%
filter(`Series ID` == sample(aus_retail$`Series ID`,1))
myseries## # A tsibble: 441 x 5 [1M]
## # Key: State, Industry [1]
## State Industry `Series ID` Month Turnover
## <chr> <chr> <chr> <mth> <dbl>
## 1 Tasmania Food retailing A3349764R 1982 Apr 28.4
## 2 Tasmania Food retailing A3349764R 1982 May 27.7
## 3 Tasmania Food retailing A3349764R 1982 Jun 27.7
## 4 Tasmania Food retailing A3349764R 1982 Jul 30.3
## 5 Tasmania Food retailing A3349764R 1982 Aug 29
## 6 Tasmania Food retailing A3349764R 1982 Sep 29.6
## 7 Tasmania Food retailing A3349764R 1982 Oct 29.5
## 8 Tasmania Food retailing A3349764R 1982 Nov 30.6
## 9 Tasmania Food retailing A3349764R 1982 Dec 35.7
## 10 Tasmania Food retailing A3349764R 1983 Jan 29.3
## # … with 431 more rows- Why is multiplicative seasonality necessary for this series?
autoplot(myseries, colour = 'blue')
Looking at the graph, we can see the seasonality of the data in the months of January, which is why multiplicative seasonality is necessary for this series.
- Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
fit <- myseries %>%
model(
`Holt-Winters’ Multiplicative` = ETS(Turnover ~ error("M") + trend("A") +
season("M")),
`Holt-Winters’ Damped Multiplicative` = ETS(Turnover ~ error("M") + trend("Ad") +
season("M"))
)
fc <- fit %>% forecast(h = "5 years")
fc %>%
autoplot(myseries, level = NULL) +
labs(title="Department Stores", subtitle = "Australian",
y="Turnover") +
guides(colour = guide_legend(title = "Forecast"))
c. Compare the RMSE of the one-step forecasts from the two methods.
Which do you prefer?
accuracy(fit) %>%select(".model", "RMSE")## # A tibble: 2 × 2
## .model RMSE
## <chr> <dbl>
## 1 Holt-Winters’ Multiplicative 3.34
## 2 Holt-Winters’ Damped Multiplicative 3.31Reviewing the previous data, we can see that the results are similar. Holt-winters Multiplicative and Holt-Winters Damped Multiplicative have very little difference. The RMSE values have a difference of 0.002167. The graphs of both forecasts are almost identical. The model with the lowest RMSE value is Holt-Winters Damped Multiplicative.
- Check that the residuals from the best method look like white noise.
fit %>% select("Holt-Winters’ Damped Multiplicative") %>% gg_tsresiduals()
Using the ACF plot, it shows that the Holt-Winters’ Damped
Multiplicative method is not white noise, because more than 5% of the
spikes are outside of the dashed lines.
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?
myseries_train <- myseries %>%
filter(year(Month) < 2011)
fit <- myseries_train %>%
model(
"Holt-Winters' Damped" = ETS(Turnover ~ error("M") + trend("Ad") +
season("M")),
"Holt-Winters' Multiplicative" = ETS(Turnover ~ error("M") + trend("A") +
season("M")),
"Seasonal Naïve Forecast" = SNAIVE(Turnover)
)
comparison <- anti_join(myseries, myseries_train,
by = c("State", "Industry", "Series ID", "Month", "Turnover"))
fc <- fit %>%
forecast(comparison)
autoplot(comparison, Turnover) +
autolayer(fc, level = NULL) +
guides(colour=guide_legend(title="Forecast")) +
ggtitle('Forecast Comparison',
subtitle= "Australian Department Stores") 
Reviewing the graph, the best model is Holt-Winters’ Damped based on the RMSE.
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?
Box-Cox transformed and find optimal lambda
lambda <- myseries_train %>%
features(Turnover, features = guerrero) %>%
pull(lambda_guerrero)
ts_bc <- myseries_train %>%
mutate(
bc_turnover = box_cox(Turnover, lambda)
)
fit <- ts_bc %>%
model(
'Box-Cox STL' = STL(bc_turnover ~ season(window = "periodic"),
robust = T),
'Box-Cox ETS' = ETS(bc_turnover)
)
best_fit <-ts_bc %>%
model(
"Holt-Winters' Damped" = ETS(Turnover ~ error("M") + trend("Ad") +
season("M"))
)
rbind(accuracy(fit),accuracy(best_fit))## # A tibble: 3 × 12
## State Indus…¹ .model .type ME RMSE MAE MPE MAPE MASE RMSSE
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Tasmania Food r… Box-C… Trai… -8.40e-4 0.0485 0.0385 -0.0226 0.641 0.300 0.322
## 2 Tasmania Food r… Box-C… Trai… -2.25e-3 0.0585 0.0474 -0.0455 0.785 0.369 0.388
## 3 Tasmania Food r… Holt-… Trai… 1.89e-1 2.93 2.22 0.165 2.54 0.383 0.392
## # … with 1 more variable: ACF1 <dbl>, and abbreviated variable name ¹IndustryAccording to the RMSE values, the Box-Cox STL model is the best performing out of the three with an RMSE of 0.04847.