HW 2 CH 8
Sal Pollastro
4/14/2022
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.
# 1 -------------------------
#View(aus_livestock)
pigs <- aus_livestock %>%
filter(Animal == "Pigs", State == "Victoria")
pigs %>%
autoplot(Count)+
labs(title="Pigs Slaughtered in Victoria") + theme(plot.title = element_text(hjust = 0.5))
# a.
ets_1a <- pigs %>%
model(ses = ETS(Count ~ error('A') + trend('N') + season('N')))
opt_values_pigs <- ets_1a %>% 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.07opt_values_pigs## # A mable: 1 x 3
## # Key: Animal, State [1]
## Animal State ses
## <fct> <fct> <model>
## 1 Pigs Victoria <ETS(A,N,N)>ets_1a## # A mable: 1 x 3
## # Key: Animal, State [1]
## Animal State ses
## <fct> <fct> <model>
## 1 Pigs Victoria <ETS(A,N,N)>pigs_f <- ets_1a %>% forecast(h = 4)
#autoplot(pigs_f) +labs(title="Four Month Forecast") +
# theme(plot.title = element_text(hjust = 0.5))
pigs_f %>%
autoplot(filter(pigs, Month >= yearmonth('2017 Jan'))) +
labs(title = 'Four Month Forecast') + theme(plot.title = element_text(hjust = 0.5))
pigs_f %>%
autoplot(filter(pigs, Month >= yearmonth('2017 Jan'))) +
labs(title = 'Four Month Forecast') + theme(plot.title = element_text(hjust = 0.5)) The optimal values are: alpha: = 0.3221247 ℓ0 = 100646.6
- 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.
# b
stdDev <- augment(ets_1a) %>% pull(.resid) %>%
sd()
yhat <- pigs_f %>%
pull(Count) %>% head(1)
upper <- yhat + 1.96 * stdDev
lower <- yhat - 1.96 * stdDev
paste(upper, lower)## [1] "N(113502, 8.7e+07) N(76871, 8.7e+07)"interval <- hilo(pigs_f$Count, 95)
interval$lower## [1] 76854.79 75927.17 75042.22 74194.54interval$upper## [1] 113518.3 114445.9 115330.9 116178.6The prediction interval is (76,871 , 113,502).
- Write your own function to implement simple exponential smoothing. The function should take arguments y (the time series), alpha (the smoothing parameter α) and level (the initial level ℓ0 ). It should return the forecast of the next observation in the series. Does it give the same forecast as ETS()?
ses2 <- function(y, alpha, level)
{
y_hat <- level
for(index in 1:length(y))
{
y_hat <- alpha*y[index] + (1 - alpha)*y_hat
}
print(paste0("My Forecast: ",
as.character(y_hat)
))
}
victorian_pigs <- aus_livestock %>%
filter(Animal == "Pigs", State == "Victoria") %>%
pull(Count)
ses_pigs_fc <- victorian_pigs %>%
ses2(alpha = 0.3221247, level = 100646.6)## [1] "My Forecast: 95186.5574351184"paste('fable :: ETS()' = pigs_f$Count[1])## [1] "N(95187, 8.7e+07)"- Modify your function from the previous exercise to return the sum of squared errors rather than the forecast of the next observation. Then use the optim() function to find the optimal values of α and ℓ0. Do you get the same values as the ETS() function?
# 3 -------------------------------------
ses2 <- function(par, y)
{
alpha <- par[1]
level <- par[2]
n <- length(y)
yhat <- numeric(n)
yhat[1] <- level
for(i in 2:n)
{
yhat[i] <- alpha * y[i-1] + (1 - alpha) * yhat[i-1]
}
return(sum((y - yhat)^2) %>%
round(3)
)
}
optParam <- optim(c(0.00001, victorian_pigs[1]), ses2, y=victorian_pigs)$par
#Optimal parameters from my SES function:
paste("Optimal parameters from my SES function:", optParam)## [1] "Optimal parameters from my SES function: 0.323696659685336"
## [2] "Optimal parameters from my SES function: 94224.3761911688"- Combine your previous two functions to produce a function that both finds the optimal values of α and ℓ0, and produces a forecast of the next observation in the series.
SES <- function(init_pars, data){
fc_next <- 0
SSE <- function(pars, data){
error <- 0
SSE <- 0
alpha <- pars[1]
l0 <- pars[2]
y_hat <- l0
for(index in 1:length(data)){
error <- data[index] - y_hat
SSE <- SSE + error^2
y_hat <- alpha*data[index] + (1 - alpha)*y_hat
}
fc_next <<- y_hat
return(SSE)
}
optim_pars <- optim(par = init_pars, data = data, fn = SSE)
return(list(
Next_observation_forecast = fc_next,
alpha = optim_pars$par[1],
l0 = optim_pars$par[2]
))
}
SES(c(0.5, victorian_pigs[1]), victorian_pigs)## $Next_observation_forecast
## [1] 95186.55
##
## $alpha
## [1] 0.3221274
##
## $l0
## [1] 99223.08- 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.
# a.
#View(global_economy)
exp_italy <- global_economy %>%
filter(Code == 'ITA')
exp_italy %>% autoplot(Exports) +
labs(title = 'Italy Exports') + theme(plot.title = element_text(hjust = 0.5))
This data series has an increasing trend throughout the years.
- Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
# b.
ets_ita <- exp_italy %>%
model(ANN = ETS(Exports ~ error("A") + trend("N") + season("N")))
ita_f <- ets_ita %>% forecast(h = 4)
ita_f %>% autoplot(exp_italy) +
labs(title = 'Italy Exports + Forecast') + theme(plot.title = element_text(hjust = 0.5))
- Compute the RMSE values for the training data.
# c.
ets_ita %>% accuracy()## # A tibble: 1 x 11
## Country .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <fct> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Italy ANN Training 0.324 1.34 1.00 1.39 4.75 0.983 0.991 -0.00701- 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.
# d.
new_ets_ita <- exp_italy %>%
model(ANN = ETS(Exports ~ error("A") + trend("N") + season("N")),
AAN = ETS(Exports ~ error("A") + trend("A") + season("N"))
)
new_ets_ita %>% accuracy() ## # A tibble: 2 x 11
## Country .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <fct> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Italy ANN Training 0.324 1.34 1.00 1.39 4.75 0.983 0.991 -0.00701
## 2 Italy AAN Training -0.00841 1.30 0.934 -0.291 4.47 0.916 0.962 0.0375The AAN model has a better RMSE value so the AAN model is a better fit.
- Compare the forecasts from both methods. Which do you think is best?
# e.
new_ets_ita_f <- new_ets_ita %>% forecast(h = 4)
new_ets_ita_f %>%
autoplot(exp_italy, level = NULL) +
labs(title = 'Italy Exports + Compared Forecasts') + theme(plot.title = element_text(hjust = 0.5))
The AAN model allows for the trend component to be included, which allows for a more accurate forecast.
- 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.
# f.
sd_5f <- new_ets_ita %>% select(Country, AAN) %>%
accuracy() %>%
transmute(Country, s = RMSE)
new_ets_ita %>%
select(Country, AAN) %>%
forecast(h = 3) %>%
left_join(sd_5f, by = 'Country') %>%
mutate(
low = Exports - 1.96 * s,
high = Exports + 1.96 * s
) %>%
select(Country, Exports, low, high)## # A fable: 3 x 5 [1Y]
## # Key: Country [1]
## Country Exports low high Year
## <fct> <dist> <dist> <dist> <dbl>
## 1 Italy N(32, 1.8) N(29, 1.8) N(34, 1.8) 2018
## 2 Italy N(32, 3.4) N(29, 3.4) N(34, 3.4) 2019
## 3 Italy N(32, 5) N(30, 5) N(35, 5) 2020- 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.
# 6 ---------------------------------------------
china_gdp <- global_economy %>%
filter(Country == "China")
china_gdp$GDP <- china_gdp$GDP/1000000000
china_gdp %>% autoplot(GDP) +
labs(title = "China GDP") +
ylab("GDP (in billions)") +
xlab("Year") + theme(plot.title = element_text(hjust = 0.5))
transformed_china_gdp <- china_gdp %>%
features(GDP, features = guerrero) %>%
pull(lambda_guerrero)
boxcox_china_gdp <- china_gdp %>% autoplot(box_cox(GDP, transformed_china_gdp))
china_ets <- china_gdp %>%
model(
ets = ETS(GDP),
ets_damped = ETS(GDP ~ trend("Ad")),
ets_boxcox = ETS(box_cox(GDP,transformed_china_gdp)),
ets_boxcox_damped = ETS(box_cox(GDP, transformed_china_gdp) ~ trend("Ad")),
ets_log = ETS(log(GDP))
)
china_ets_f <- china_ets %>% forecast(h = "20 years")
china_ets_f %>% autoplot(china_gdp, level = NULL) +
labs(title = "Forecasted China GDP") +
ylab("GDP (in billions)") +
xlab("Year") + theme(plot.title = element_text(hjust = 0.5))
The ETS and ETS boxcox forecasts seem like the upper end of a forecast and would be unlikely. The ETS boxcox damped seems like it shows an increasing exponential trend. The ETS Damped and ETS Log seem like they are more accurate forecasts.
- 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?
# 7 ---------------------------------------
#View(aus_production)
aus_production %>% autoplot(Gas) +
labs(title = "Australian Gas Production") +
ylab("Production (in bcf)") +
xlab("Quarter") + theme(plot.title = element_text(hjust = 0.5))
gas_ts <- ts(aus_production$Gas, frequency = 4)
gas_ets <- ets(gas_ts, model = "MAM")
gas_f <- gas_ets %>% forecast(h=12)
gas_f %>% autoplot() +
labs(title = "Australian Gas Production with Forecast") +
ylab("Production (in bcf)") +
xlab("Quarter") + theme(plot.title = element_text(hjust = 0.5))
gas_ets_damped <- ets(gas_ts, model = "MAM", damped = TRUE)
gas_ets_damped_f <- gas_ets_damped %>% forecast(h = 12)
gas_ets_damped_f %>% autoplot() +
labs(title = "Australian Gas Production with Damped Forecast") +
ylab("Production (in bcf)") +
xlab("Quarter") + theme(plot.title = element_text(hjust = 0.5))
gas_ets## ETS(M,A,M)
##
## Call:
## ets(y = gas_ts, model = "MAM")
##
## Smoothing parameters:
## alpha = 0.6529
## beta = 0.1442
## gamma = 0.0978
##
## Initial states:
## l = 5.9456
## b = 0.0706
## s = 0.9309 1.1779 1.0749 0.8163
##
## sigma: 0.057
##
## AIC AICc BIC
## 1680.929 1681.794 1711.389gas_ets_damped## ETS(M,Ad,M)
##
## Call:
## ets(y = gas_ts, model = "MAM", damped = TRUE)
##
## Smoothing parameters:
## alpha = 0.6489
## beta = 0.1551
## gamma = 0.0937
## phi = 0.98
##
## Initial states:
## l = 5.8589
## b = 0.0994
## s = 0.9282 1.1779 1.0768 0.8171
##
## sigma: 0.0573
##
## AIC AICc BIC
## 1684.028 1685.091 1717.873Multiplicative seasonality is necessary because of the increasing seasonal variation.
- Recall your retail time series data (from Exercise 8 in Section 2.10).
- Why is multiplicative seasonality necessary for this series?
# a.
set.seed(12345678)
myseries <- aus_retail %>%
filter(`Series ID` == sample(aus_retail$`Series ID`,1))
retail_ts <- ts(myseries$Turnover, frequency = 12)
retail_ts %>% autoplot() +
labs(title = "Australian Retail Turnover") +
ylab("Turnover") +
theme(plot.title = element_text(hjust = 0.5))
Multiplicative seasonality is necessary because of the increasing seasonal variation.
- Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
# b.
retail_holtwinters <- ets(retail_ts, model = "MAM")
retail_f <- retail_holtwinters %>% forecast(h=36)
retail_f %>% autoplot() +
labs(title = "Australian Retail Turnover with Forecast") +
ylab("Turnover") +
theme(plot.title = element_text(hjust = 0.5))
retail_holtwinters_damped <- ets(retail_ts, model = "MAM", damped = TRUE)
retail_hw_damped_f <- retail_holtwinters_damped %>% forecast(h=36)
retail_hw_damped_f %>% autoplot() +
labs(title = "Australian Retail Turnover with Damped Forecast") +
ylab("Turnover") +
theme(plot.title = element_text(hjust = 0.5))
The damped forecasts produce flatter forecasts that probably don’t account for the trend component.
- Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
# c.
sqrt(mean(retail_f$residuals^2))## [1] 0.06537331sqrt(mean(retail_hw_damped_f$residuals^2))## [1] 0.06600816- Check that the residuals from the best method look like white noise.
# d.
checkresiduals(retail_f)
##
## Ljung-Box test
##
## data: Residuals from ETS(M,A,M)
## Q* = 34.479, df = 8, p-value = 3.326e-05
##
## Model df: 16. Total lags used: 24- 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)
myseriesfit <- 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")),
"SNAIVE Forecast" = SNAIVE(Turnover)
)
Compare <- anti_join(myseries, myseries_train,
by = c("State", "Industry", "Series ID", "Month", "Turnover"))
fc <- myseriesfit %>%
forecast(Compare)
autoplot(Compare, Turnover) +
autolayer(fc, level = NULL) +
guides(colour=guide_legend(title="Forecast")) +
labs(title = "Forecast Comparison") +
theme(plot.title = element_text(hjust = 0.5))
- 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?
# 9 ---------------------------------------------
lambda <- myseries_train %>%
features(Turnover, features = guerrero) %>%
pull(lambda_guerrero)
ts_bc <- myseries_train %>%
mutate(
bc = box_cox(Turnover, lambda)
)
bc_fit <- ts_bc %>%
model(
'STL (BoxCox)' = STL(bc ~ season(window = "periodic"),
robust = T),
'ETS (BoxCox)' = ETS(bc)
)
bc_fit_opt <-ts_bc %>%
model(
"Holt-Winters' Multiplicative" = ETS(Turnover ~ error("M") +
trend("A") +
season("M"))
)
rbind(accuracy(bc_fit),accuracy(bc_fit_opt))## # A tibble: 3 x 12
## State Industry .model .type ME RMSE MAE MPE MAPE MASE RMSSE
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 North~ Clothing,~ STL (~ Trai~ 0.00530 0.0687 0.0524 0.125 2.63 0.325 0.319
## 2 North~ Clothing,~ ETS (~ Trai~ 0.00668 0.0832 0.0654 0.258 3.25 0.405 0.387
## 3 North~ Clothing,~ Holt-~ Trai~ -0.0119 0.518 0.384 -0.446 5.21 0.420 0.427
## # ... with 1 more variable: ACF1 <dbl>- Compute the total domestic overnight trips across Australia from the tourism dataset.
- Plot the data and describe the main features of the series.
# a.
aus_trip <- tourism %>%
summarise(Trips = sum(Trips))
aus_trip %>% autoplot(Trips) +
labs(title = "Australian Overnight Trips") +
theme(plot.title = element_text(hjust = 0.5))
- Decompose the series using STL and obtain the seasonally adjusted data.
# b.
aus_trip_stl <- aus_trip %>%
model(STL(Trips)) %>% components()
aus_trip_stl %>% autoplot()
aus_trip_stl %>% as_tsibble() %>%
autoplot(season_adjust) +
labs(title = "Australian Overnight Trips (STL)") +
theme(plot.title = element_text(hjust = 0.5))
- Forecast the next two years of the series using an additive damped trend method applied to the seasonally adjusted data. (This can be specified using decomposition_model().)
# c.
aus_trip %>%
model(
decomposition_model(STL(Trips),
ETS(season_adjust ~ error("A") + trend("Ad") + season("N")
)
)
) %>%
forecast(h = "2 years") %>% autoplot(aus_trip) +
labs(title = "Australian Overnight Trips with Additive Damped Trend Forecast") +
theme(plot.title = element_text(hjust = 0.5))
- Forecast the next two years of the series using an appropriate model for Holt’s linear method applied to the seasonally adjusted data (as before but without damped trend).
# d.
aus_trip %>%
model(
decomposition_model(STL(Trips),
ETS(season_adjust ~ error("A") +
trend("A") +
season("N")
)
)
) %>%
forecast(h = "2 years") %>%
autoplot(aus_trip) +
labs(title = "Australian Overnight Trips with Holt's Linear Method Forecast") +
theme(plot.title = element_text(hjust = 0.5))
- Now use ETS() to choose a seasonal model for the data.
# e.
aus_trip %>%
model(
ETS(Trips)
) %>%
forecast(h = "2 years") %>%
autoplot(aus_trip) +
labs(title = "Australian Overnight Trips (ETS)") +
theme(plot.title = element_text(hjust = 0.5))
- Compare the RMSE of the ETS model with the RMSE of the models you obtained using STL decompositions. Which gives the better in-sample fits?
# f.
compare_STL <- aus_trip %>%
model(
STL_AAdN = decomposition_model(STL(Trips),
ETS(season_adjust ~ error("A") +
trend("Ad") +
season("N"))),
STL_AAN = decomposition_model(STL(Trips),
ETS(season_adjust ~ error("A") +
trend("A") +
season("N"))),
ETS = ETS(Trips)
)
accuracy(compare_STL)## # A tibble: 3 x 10
## .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 STL_AAdN Training 103. 763. 576. 0.367 2.72 0.607 0.629 -0.0174
## 2 STL_AAN Training 99.7 763. 574. 0.359 2.71 0.604 0.628 -0.0182
## 3 ETS Training 105. 794. 604. 0.379 2.86 0.636 0.653 -0.00151The STL AAN model performed the best based on the lower diagnostic statistics.
- Compare the forecasts from the three approaches? Which seems most reasonable?
# g.
compare_STL %>%
forecast(h = "2 years") %>%
autoplot(aus_trip, level = NULL) +
labs(title = "Australian Overnight Trips: Comparing Forecasts") +
guides(colour=guide_legend(title="Types of Forecast")) +
theme(plot.title = element_text(hjust = 0.5))
The model with the lowest Root Squared Estimate value would be the STL AAN model.
- Check the residuals of your preferred model.
# h.
resid_aan <- compare_STL %>% select(STL_AAN)
resid_aan %>% gg_tsresiduals()## Warning: Removed 4 row(s) containing missing values (geom_path).## Warning: Removed 4 rows containing missing values (geom_point).## Warning: Removed 4 rows containing non-finite values (stat_bin).
- For this exercise use the quarterly number of arrivals to Australia from New Zealand, 1981 Q1 – 2012 Q3, from data set aus_arrivals.
- Make a time plot of your data and describe the main features of the series.
# a.
#View(aus_arrivals)
aus_nz_arrivals <- aus_arrivals %>% filter(Origin == 'NZ')
autoplot(aus_nz_arrivals, Arrivals) +
labs(title = "New Zealand Arrivals in Australia") +
theme(plot.title = element_text(hjust = 0.5))
This time series has an increasing trend and varying seasonality over the years.
- Create a training set that withholds the last two years of available data. Forecast the test set using an appropriate model for Holt-Winters’ multiplicative method.
# b.
arrivals_train <- aus_nz_arrivals %>%
slice(1:(nrow(aus_nz_arrivals)-8))
arrivals_train %>%
model(
ETS(Arrivals ~ error("M") +
trend("A") +
season("M"))
) %>%
forecast(h = "2 years") %>%
autoplot(level = NULL, colour = "red") +
autolayer(aus_nz_arrivals, Arrivals) +
labs(title = "New Zealand Arrivals in Australia: Train vs. Actual") +
theme(plot.title = element_text(hjust = 0.5))
Why is multiplicative seasonality necessary here? Multiplicative seasonality is necessary because of the increasing seasonal variation.
Forecast the two-year test set using each of the following methods: an ETS model; an additive ETS model applied to a log transformed series; a seasonal naïve method; an STL decomposition applied to the log transformed data followed by an ETS model applied to the seasonally adjusted (transformed) data.
# d.
compare_models <- anti_join(aus_nz_arrivals,
arrivals_train,
by = c("Quarter",
"Origin",
"Arrivals"))
arrivals_models <- aus_nz_arrivals %>%
model(
ETS_mod = ETS(Arrivals),
Log_mod = ETS(log(Arrivals) ~ error("A") +
trend("A") +
season("A")),
SNAIVE_mod = SNAIVE(Arrivals),
STL_decomp_log = decomposition_model(STL(log(Arrivals)),
ETS(season_adjust))
)
arrivals_f <- arrivals_models %>% forecast(h = "2 years")
arrivals_f %>% autoplot(level = NULL) +
autolayer(compare_models, Arrivals) +
labs(title = "New Zealand Arrivals in Australia: Train vs. Actual") +
theme(plot.title = element_text(hjust = 0.5)) +
guides(colour = guide_legend(title = "Forecast"))
- Which method gives the best forecasts? Does it pass the residual tests?
# e.
aus_optimal <- arrivals_models %>% select('ETS_mod')
aus_optimal %>% gg_tsresiduals()
augment(aus_optimal) %>% features(.resid, ljung_box)## # A tibble: 1 x 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 ETS_mod 0.596 0.440the ETS model is the best forecase as the lags in the ACF plot lie within the bounds.
- Compare the same four methods using time series cross-validation instead of using a training and test set. Do you come to the same conclusions?
# f.
arrivals_models_crossval <- aus_nz_arrivals %>%
slice(1:(n()-3)) %>%
stretch_tsibble(.init = 36, .step = 3)
arrivals_models_crossval %>%
model(
ETS_mod = ETS(Arrivals),
Add_Log = ETS(log(Arrivals) ~ error("A") +
trend("A") +
season("A")),
SNAIVE_mod = SNAIVE(Arrivals),
STL_decomp_log = decomposition_model(STL(log(Arrivals)),
ETS(season_adjust))
) %>%
forecast(h = 3) %>%
accuracy(aus_nz_arrivals)## # A tibble: 4 x 11
## .model Origin .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Add_Log NZ Test -375. 14347. 11021. 0.200 5.83 0.741 0.746 0.309
## 2 ETS_mod NZ Test 4627. 15327. 11799. 2.23 6.45 0.793 0.797 0.283
## 3 SNAIVE_mod NZ Test 8244. 18768. 14422. 3.83 7.76 0.970 0.976 0.566
## 4 STL_decomp_log NZ Test 4252. 15618. 11873. 2.04 6.25 0.798 0.812 0.244The additive log model is the best model, not the ETS model as before, because the RSME value is lower this time for the additive log.
- Apply cross-validation techniques to produce 1 year ahead ETS and seasonal naïve forecasts for Portland cement production (from aus_production). Use a stretching data window with initial size of 5 years, and increment the window by one observation.
# a.
#View(aus_production)
portland_cement_crossval <- aus_production %>%
slice(1:(n()-4)) %>%
stretch_tsibble(.init = 20, .step = 1)
port_cement_crossval_f <- portland_cement_crossval %>%
model(ETS(Cement),
SNAIVE(Cement)
) %>%
forecast(h = "1 year")- Compute the MSE of the resulting 4-step-ahead errors. Comment on which forecasts are more accurate. Is this what you expected?
port_cement_crossval_f %>%
group_by(.id, .model) %>%
mutate(h = row_number()) %>%
ungroup() %>%
accuracy(aus_production, by = c(".model", "h"))## # A tibble: 8 x 11
## .model h .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <int> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ETS(Cement) 1 Test -0.0902 82.9 60.1 -0.227 3.95 0.597 0.625 -0.00185
## 2 ETS(Cement) 2 Test -0.653 101. 72.0 -0.325 4.74 0.708 0.756 0.495
## 3 ETS(Cement) 3 Test -1.71 119. 87.0 -0.492 5.80 0.856 0.894 0.616
## 4 ETS(Cement) 4 Test -0.729 137. 102. -0.543 6.65 1.01 1.03 0.699
## 5 SNAIVE(Ceme~ 1 Test 30.9 138. 107. 1.97 6.99 1.06 1.04 0.640
## 6 SNAIVE(Ceme~ 2 Test 30.0 139. 107. 1.90 6.96 1.05 1.04 0.649
## 7 SNAIVE(Ceme~ 3 Test 29.5 139. 107. 1.85 6.95 1.05 1.04 0.651
## 8 SNAIVE(Ceme~ 4 Test 30.8 140. 108 1.91 6.99 1.06 1.05 0.637- Compare ETS(), SNAIVE() and decomposition_model(STL, ???) on the following five time series. You might need to use a Box-Cox transformation for the STL decomposition forecasts. Use a test set of three years to decide what gives the best forecasts.
Beer and bricks production from aus_production.
aus_production %>% autoplot(Beer) +
labs(title = "Australian Beer Production") +
theme(plot.title = element_text(hjust = 0.5))
beer_f <- aus_production %>%
slice(1:(nrow(aus_production)-12)) %>%
model(
ETS_beer = ETS(Beer),
SNAIVE_beer = SNAIVE(Beer),
STL_Decomp_beer = decomposition_model(STL(log(Beer)),
ETS(season_adjust))
) %>% forecast(h = "3 years")
beer_f %>%
autoplot(filter_index(aus_production, "2000 Q1" ~ .), level = NULL) +
labs(title = "Australian Beer Production") +
theme(plot.title = element_text(hjust = 0.5))
beer_f %>% accuracy(aus_production)## # A tibble: 3 x 10
## .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ETS_beer Test 0.127 9.62 8.92 0.00998 2.13 0.563 0.489 0.376
## 2 SNAIVE_beer Test -2.92 10.8 9.75 -0.651 2.34 0.616 0.549 0.325
## 3 STL_Decomp_beer Test -2.85 9.87 8.95 -0.719 2.16 0.565 0.502 0.283Cost of drug subsidies for diabetes (ATC2 == “A10”) and corticosteroids (ATC2 == “H02”) from PBS.
# Bricks
aus_bricks <- aus_production %>% filter(!is.na(Bricks)) #### FIXXXXXXX
aus_bricks %>% autoplot(Bricks) +
labs(title = "Australian Quarterly Brick Production") +
theme(plot.title = element_text(hjust = 0.5))
bricks_f <- aus_bricks %>% slice(1:(nrow(aus_bricks)-12)) %>%
model(
ETS = ETS(Bricks),
SNAIVE = SNAIVE(Bricks),
STL_Decomp = decomposition_model(STL(log(Bricks)),
ETS(season_adjust))
) %>%
forecast(h = "3 years")
bricks_f %>%
autoplot(filter_index(aus_bricks, "1980 Q1" ~ .), level = NULL) +
labs(title = "Australian Quarterly Brick Production") +
theme(plot.title = element_text(hjust = 0.5))
bricks_f %>% accuracy(aus_bricks)## # A tibble: 3 x 10
## .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ETS Test 2.27 17.5 13.2 0.474 3.31 0.365 0.354 0.339
## 2 SNAIVE Test 32.6 36.5 32.6 7.85 7.85 0.898 0.739 -0.322
## 3 STL_Decomp Test 9.71 18.7 14.9 2.29 3.65 0.411 0.378 0.0933Total food retailing turnover for Australia from aus_retail.
# Drug Subsidies
#View(PBS)
drug_subsidies <- PBS %>%
filter(ATC2 %in% c("A10", "H02")) %>%
group_by(ATC2) %>%
summarise(Cost = sum(Cost))
drug_subsidies %>%
autoplot(Cost) + facet_wrap(vars(ATC2, scales = "free_y")) +
labs(title = "Diabetes & Corticosteroids Drug Subsidy Costs") +
theme(plot.title = element_text(hjust = 0.5))
#library(latticeExtra)
#doubleYScale(drug_subsidies$ATC2,drug_subsidies$Cost, add.ylab2 = TRUE, use.style = FALSE)
cortico <- drug_subsidies %>%
filter(ATC2 %in% 'H02')
cortico_cost <- cortico %>%
features(Cost, features = guerrero) %>%
pull(lambda_guerrero)
cortico_f <-cortico %>% filter(Month < max(Month) - 35) %>%
model(
ETS = ETS(Cost),
SNAIVE = SNAIVE(Cost),
STL_Decomp = decomposition_model(STL(box_cox(log(Cost), cortico_cost)), ETS(season_adjust))
) %>%
forecast(h = "3 years")
diabetes <- drug_subsidies %>%
filter(ATC2 %in% 'A10')
diabetes_cost <- diabetes %>%
features(Cost, features = guerrero) %>%
pull(lambda_guerrero)
diabetes_f <- diabetes %>%
filter(Month < max(Month) - 35) %>%
model(
ETS = ETS(Cost),
SNAIVE = SNAIVE(Cost),
STL_Decomp = decomposition_model(STL(box_cox(log(Cost), diabetes_cost)),
ETS(season_adjust))
) %>%
forecast(h = "3 years")
compare_drug_costs <- bind_rows(cortico_f, diabetes_f)
compare_drug_costs %>%
autoplot(drug_subsidies, level = NULL) +
labs(title = "Diabetes & Corticosteroids Drug Subsidy Costs with Forecasts") +
theme(plot.title = element_text(hjust = 0.5))
#View(aus_retail)
aus_food <- aus_retail %>%
filter(Industry=="Food retailing") %>%
summarise(Turnover=sum(Turnover))
aus_food %>%
autoplot(Turnover) +
labs(title = "Australian Food Turnover") +
theme(plot.title = element_text(hjust = 0.5))
food_turnover <- aus_food %>% features(Turnover, features = guerrero) %>%
pull(lambda_guerrero)
aus_food_f <- aus_food %>%
filter(Month < max(Month) - 35) %>%
model(ETS = ETS(Turnover),
SNAIVE = SNAIVE(Turnover),
STL_Decomp = decomposition_model(STL(box_cox(log(Turnover),
food_turnover)),
ETS(season_adjust))
) %>%
forecast(h = "3 years")
aus_food_f %>% autoplot(filter_index(aus_food, "2005 Jan" ~ .), level = NULL) +
labs(title = "Australian Food Turnover Forecast Comparisons") +
theme(plot.title = element_text(hjust = 0.5))
aus_food_f %>% accuracy(aus_food)## # A tibble: 3 x 10
## .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ETS Test -151. 194. 170. -1.47 1.65 0.639 0.634 0.109
## 2 SNAIVE Test 625. 699. 625. 5.86 5.86 2.35 2.29 0.736
## 3 STL_Decomp Test -8.34 155. 132. -0.135 1.24 0.496 0.506 0.256- Use ETS() to select an appropriate model for the following series: total number of trips across Australia using tourism, the closing prices for the four stocks in gafa_stock, and the lynx series in pelt. Does it always give good forecasts?
# a.
aus_tourism <- tourism %>%
summarise(Trips = sum(Trips))
tourism_trips <- aus_tourism %>%
model(ETS(Trips))
tourism_trips %>%
forecast() %>%
autoplot(aus_tourism) +
labs(title = "Australian Tourism Forecast") +
theme(plot.title = element_text(hjust = 0.5))
gafa_stock %>%
autoplot(Close) +
facet_grid(vars(Symbol), scales = "free_y") +
labs(title = "Closing Stock Prices") +
theme(plot.title = element_text(hjust = 0.5))
gafa_avg <- gafa_stock %>%
group_by(Symbol) %>%
mutate(trading_day = row_number()) %>%
ungroup() %>%
as_tsibble(index = trading_day, regular = TRUE)
gafa_fc <- fit <- gafa_avg %>%
model(
ETS(Close)
)
gafa_fc %>%
forecast(h=50) %>%
autoplot(gafa_avg %>% group_by_key() %>% slice((n() - 100):n())) +
labs(title = "Forecasted Closing Stock Prices") +
theme(plot.title = element_text(hjust = 0.5))## `mutate_if()` ignored the following grouping variables:
## Column `Symbol`
pelt %>%
model(ETS(Lynx))## # A mable: 1 x 1
## `ETS(Lynx)`
## <model>
## 1 <ETS(A,N,N)>pelt %>%
model(ETS(Lynx)) %>% forecast(h = 10) %>%
autoplot(pelt) +
labs(title = "Forecasted Pelt Trades") +
theme(plot.title = element_text(hjust = 0.5))
The ETS does not always give good forecasts as seen by the pelt forecast, where the forecast does not seem to give any accurate indication of future pelt trades.
- Show that the point forecasts from an ETS(M,A,M) model are the same as those obtained using Holt-Winters’ multiplicative method.
# 15 --------------------------------------------------------
usdeaths.ets2 <- ets(usdeaths,model = "MAM")
usdeaths.ets2F <- forecast(usdeaths.ets2, h=1)
usdeaths.hwM <- hw(usdeaths, seasonal = 'multiplicative', h=1)
usdeaths.ets2F$mean## Jan
## 1979 8233.107usdeaths.hwM$mean## Jan
## 1979 8217.64- Show that the forecast variance for an ETS(A,N,N) model is given by σ^2[1+α^2(h−1)]
ibm_ets <- ets(ibmclose, model = "ANN")
ibm_ets_fc <- forecast(ibm_ets, 1)
ibm_ets_fc## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 370 356.9995 347.6907 366.3083 342.7629 371.2361ibm_sigma <- 7.2637
alpha <- .9999^2
h <- 2
conf_int <- ibm_sigma*(1+alpha*(h-1))
ibm_ets_fc$mean[1]## [1] 356.9995ibm_ets_fc$mean[1] - conf_int## [1] 342.4736ibm_ets_fc$lower[1, '95%']## 95%
## 342.7629- Write down 95% prediction intervals for an ETS(A,N,N) model as a function of ℓT, α, h and σ, assuming normally distributed errors.
^yT+h|T ± 1.96^estimate_std_dev