HW5_624 Predictive Analytics
Karim Hammoud
8.1
Consider the the number of pigs slaughtered in Victoria, available in the aus_livestock dataset.
head(aus_livestock)## # A tsibble: 6 x 4 [1M]
## # Key: Animal, State [1]
## Month Animal State Count
## <mth> <fct> <fct> <dbl>
## 1 1976 Jul Bulls, bullocks and steers Australian Capital Territory 2300
## 2 1976 Aug Bulls, bullocks and steers Australian Capital Territory 2100
## 3 1976 Sep Bulls, bullocks and steers Australian Capital Territory 2100
## 4 1976 Oct Bulls, bullocks and steers Australian Capital Territory 1900
## 5 1976 Nov Bulls, bullocks and steers Australian Capital Territory 2100
## 6 1976 Dec Bulls, bullocks and steers Australian Capital Territory 1800x <- aus_livestock %>%
filter(Animal=="Pigs", State=="Victoria")
x %>% autoplot(Count)
- 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 <- x %>% model(ses=ETS(Count~error("A")+trend("N") +season("N")))
report <- fit %>% report()## 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.07We can see a = 0.3221247 and ℓ0 = 100646.6
Now I will generate forecasts for the next four months
forc <- fit %>% forecast(h="4 months")
forc## # 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 ses 2019 Jan N(95187, 8.7e+07) 95187.
## 2 Pigs Victoria ses 2019 Feb N(95187, 9.7e+07) 95187.
## 3 Pigs Victoria ses 2019 Mar N(95187, 1.1e+08) 95187.
## 4 Pigs Victoria ses 2019 Apr N(95187, 1.1e+08) 95187.For the above table we can see the forcast for the next four months.
Lets plots with the last 8 months.
forc %>%
autoplot(filter(aus_livestock, Month >= yearmonth("2018 May"))) +
labs(title="forecast")
- 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.
s <- augment(fit) %>% pull(.resid) %>%
sd()
yhat <- forc %>%
pull(Count) %>%
head(1)
yhat + c(-1,1) * 1.96 *s## <distribution[2]>
## [1] N(76871, 8.7e+07) N(113502, 8.7e+07)The prediction interval is 76871≤ŷ≤113502.
Now I will generate the prediction interval produced by R
forc %>%
mutate(interval = hilo(Count, 95)) %>%
pull(interval)## <hilo[4]>
## [1] [76854.79, 113518.3]95 [75927.17, 114445.9]95 [75042.22, 115330.9]95
## [4] [74194.54, 116178.6]95The itnerval look close to each other, because the R cacluated teh variance of the residuals in a different way.
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.
BEL <- global_economy %>% filter(Code=="BEL")
BEL %>% autoplot(Exports)
It looks the exports has a upward trend.
- Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
fit <- BEL %>%
model(AMM = ETS(Exports~ error("A") + trend("N") + season ("N")))
forc <- fit %>% forecast(h=4)
forc %>% autoplot(BEL)
- Compute the RMSE values for the training data.
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 Belgium AMM Training 0.877 3.23 2.52 1.35 4.30 0.987 0.989 -0.0642The RMSE value is 3.227668.
- 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.
fit_compare <- BEL %>%
model(ANN = ETS(Exports ~ error("A") + trend ("N") + season("N")),
AAN = ETS(Exports ~ error("A") + trend ("N") + season("N")))
accuracy(fit_compare)## # 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 Belgium ANN Training 0.877 3.23 2.52 1.35 4.30 0.987 0.989 -0.0642
## 2 Belgium AAN Training 0.877 3.23 2.52 1.35 4.30 0.987 0.989 -0.0642AAN reduces the RMSE value, which results into better accuracy.
- Compare the forecasts from both methods. Which do you think is best?
fit_compare %>%
forecast(h=4) %>%
autoplot(BEL, level = NULL)
AAN model is higher than ANN model because AAN model allows increasing trend, so AAN model is the best model for trend models.
- 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.
Using AAN because its better I think.
s <- fit_compare %>%
select(Country, AAN) %>%
accuracy() %>%
transmute(Country, s = RMSE)
fit_compare %>%
select(Country, AAN) %>%
forecast(h=1) %>%
left_join(s, by = "Country") %>%
mutate( low = Exports - 1.96 * s,
high = Exports + 1.96 * s) %>%
select(Country, Exports, low, high)## # A fable: 1 x 5 [?]
## # Key: Country [1]
## Country Exports low high Year
## <fct> <dist> <dist> <dist> <dbl>
## 1 Belgium N(85, 11) N(79, 11) N(91, 11) 20188.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)
Because of the increase in the vairance, I will use BOX COX transformation to solve this problem.
lambda <- china %>%
features(GDP, features = guerrero) %>%
pull(lambda_guerrero)
bc_china <- china %>% autoplot(box_cox(GDP, lambda))
head(bc_china)## $data
## # A tsibble: 58 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 China CHN 1960 59716467625. NA NA 4.43 4.31 667070000
## 2 China CHN 1961 50056868958. -27.3 NA 3.49 3.87 660330000
## 3 China CHN 1962 47209359006. -5.58 NA 2.91 4.05 665770000
## 4 China CHN 1963 50706799903. 10.3 NA 2.86 4.01 682335000
## 5 China CHN 1964 59708343489. 18.2 NA 2.86 3.77 698355000
## 6 China CHN 1965 70436266147. 17.0 NA 3.19 3.64 715185000
## 7 China CHN 1966 76720285970. 10.7 NA 3.24 3.49 735400000
## 8 China CHN 1967 72881631327. -5.77 NA 2.98 3.28 754550000
## 9 China CHN 1968 70846535056. -4.10 NA 2.92 3.30 774510000
## 10 China CHN 1969 79705906247. 16.9 NA 2.41 3.05 796025000
## # … with 48 more rows
##
## $layers
## $layers[[1]]
## geom_line: na.rm = FALSE, orientation = NA
## stat_identity: na.rm = FALSE
## position_identity
##
##
## $scales
## <ggproto object: Class ScalesList, gg>
## add: function
## clone: function
## find: function
## get_scales: function
## has_scale: function
## input: function
## n: function
## non_position_scales: function
## scales: list
## super: <ggproto object: Class ScalesList, gg>
##
## $mapping
## Aesthetic mapping:
## * `x` -> `Year`
## * `y` -> `box_cox(GDP, lambda)`
##
## $theme
## list()
##
## $coordinates
## <ggproto object: Class CoordCartesian, Coord, gg>
## aspect: function
## backtransform_range: function
## clip: on
## default: TRUE
## distance: function
## expand: TRUE
## is_free: function
## is_linear: function
## labels: function
## limits: list
## modify_scales: function
## range: function
## render_axis_h: function
## render_axis_v: function
## render_bg: function
## render_fg: function
## setup_data: function
## setup_layout: function
## setup_panel_guides: function
## setup_panel_params: function
## setup_params: function
## train_panel_guides: function
## transform: function
## super: <ggproto object: Class CoordCartesian, Coord, gg>Not it looks more of a trend model.
fit <- china %>%
model(est=ETS(GDP), #ETS model
ets_damped = ETS(GDP ~ trend("ad")), #damped model
ets_bc = ETS(box_cox(GDP, lambda)), # ets with box cox trasnformation
ets_bc_damped = ETS(box_cox(GDP, lambda) ~ trend ("Ad")) ,#damped with boxcox trasnforamtion
ets_log = ETS(log(GDP))
)
fit %>%
forecast(h="10 years") %>%
autoplot(china, level = NULL)
For the chart we can see the box-cox and log have great impact on the forecast, while damping has a small effect.
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)
fit <- aus_production %>%
filter(Quarter > yearquarter("1990 Q4")) %>%
model(fit = ETS(Gas))
report(fit)## Series: Gas
## Model: ETS(A,A,A)
## Smoothing parameters:
## alpha = 0.4902673
## beta = 0.0001000091
## gamma = 0.0001001453
##
## Initial states:
## l[0] b[0] s[0] s[-1] s[-2] s[-3]
## 133.4487 1.227524 -12.80569 26.37931 10.65398 -24.2276
##
## sigma^2: 28.4997
##
## AIC AICc BIC
## 610.6743 613.3213 631.8847fit %>%
forecast(h = "3 years") %>%
autoplot(aus_production)
The multiplicative seasonality is necessary because the seasonal variation increases overtime.
8.8
Recall your retail time series data (from Exercise 8 in Section 2.10).
set.seed(123)
myseries <- aus_retail %>%
filter(`Series ID` == "A3349337W")
myseries %>% autoplot(Turnover)
- Why is multiplicative seasonality necessary for this series?
Because for the chart we can see there is an increase over time, so multiplicative seasonality necessary is really important because of the increase is not consistant.
- Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
fit <- myseries %>% model( "Holt-Winters' Damped" = ETS(Turnover ~ error("M") + trend("Ad") + season("M")),
"Holt-Winters' Multiplicative" = ETS(Turnover ~ error("M") + trend("A") + season("M")))
fc <- fit %>% forecast(h = "10 years")
fc %>%
autoplot(myseries, level = NULL) + labs(title="Australian retail", y="Turnover") +
guides(colour = guide_legend(title = "Forecast"))
The damped methoed is less than the mutiplicative.
- 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' Damped 14.6
## 2 Holt-Winters' Multiplicative 15.0When comparing the two methods the that Holt-Winters’ damped model has better RMSE than Holt-Winters’ multiplicative method. So Holt-Winter’s damped is better for this time series regarding the RMSE value.
- Check that the residuals from the best method look like white noise.
fit %>% select("Holt-Winters' Damped") %>% gg_tsresiduals()
THE plot and histogram, we can say the residuals look like a white noise with occasional spikes, but when we check the ACF plot, it shows this method is not white noise, because more than 5% of the spikes is out of the bound. So we will try another method (multiplicative)
fit %>% select("Holt-Winters' Multiplicative") %>% gg_tsresiduals()
Multiplicative yield plots are better than damped. We can see on the ACF plot, spikes that out of the bound less than before.
- 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?
# Split the data we use
myseries_train <- myseries %>%
filter(year(Month) < 2011)
# Add seasonal naïve model to our fit model
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))## Warning: 1 error encountered for Holt-Winters' Damped
## [1] Invalid error type# Assign the data over 2011 to a variable as our forecast comparison later
Compare <- anti_join(myseries, myseries_train, by = c("State", "Industry", "Series ID", "Month", "Turnover"))
# Do the forecasting according to comparison data
fc <- fit %>% forecast(Compare)
# Call the Comparison plot
autoplot(Compare, Turnover) +
autolayer(fc, level = NULL) + guides(colour=guide_legend(title="Forecast")) + ggtitle('Comparison of Forecast')## Warning: Removed 96 row(s) containing missing values (geom_path).
We see Holt-Winters’ Multiplicative method can beat the seasonal naïve approach. Otherwise Holt- Winters’Damped method didn’t show any significant difference with seasonal naïve approach.
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 <- myseries_train %>%
features(Turnover, features = guerrero) %>%
pull(lambda_guerrero)
ts_bc <- myseries_train %>%
mutate(bc = box_cox(Turnover, lambda))
fit <- ts_bc %>%
model("STL (BoxCox)" = STL(bc~season(window = "periodic"),
robust = T),
'ETS (BoxCox)' = ETS(bc))
fit_best <- ts_bc %>%
model("Holt-Winters' Multiplicative" = ETS(Turnover ~ error("M")+ trend("A")+ season("M"))
)
rbind(accuracy(fit), accuracy(fit_best))## # A tibble: 3 × 12
## State Industry .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 New S… Hardwar… STL (… Trai… 0.320 12.9 8.67 -0.399 5.40 0.421 0.507 0.355
## 2 New S… Hardwar… ETS (… Trai… -1.74 13.2 9.37 -1.47 5.85 0.455 0.519 0.221
## 3 New S… Hardwar… Holt-… Trai… -1.14 13.5 9.55 -0.912 5.79 0.450 0.515 0.247We can see the box cox trasnformation is better