Homework 5
Nikoleta Emanouilidi
2024-10-03
library(fpp3)
library(fable)
library(ggplot2)
library(tsibble)
library(dplyr)Question 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 a and l0 , and generate forecasts for the next four months.
pigs_victoria <- aus_livestock |>
filter(Animal=='Pigs', State=='Victoria')|>
model(ETS(Count ~ error("A") + trend("N") + season("N")))
report(pigs_victoria)## 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.07forec <- pigs_victoria |>
forecast(h = 4)
forec |>
autoplot(aus_livestock) +
labs(y="Count", title="Number of pigs slaughtered in Victoria") +
guides(colour = "none")
I found that the optimal values of alpha is 0.3221247 and for l[0] is 100646.6
b) Compute a 95% prediction interval for the first forecast using yhat +- 1.96s, where s is the standard deviation of the residuals. Compare your interval with the interval produced by R.
## Generated by R
forec %>% hilo() %>% pull(-1) ## <hilo[4]>
## [1] [76854.79, 113518.3]95 [75927.17, 114445.9]95 [75042.22, 115330.9]95
## [4] [74194.54, 116178.6]95## Calculated by hand
residuals_std <- sd(augment(pigs_victoria)$`.resid`)
# Compute confidence interval
lower_bound <- forec$.mean[1] - (residuals_std * 1.96)
upper_bound <- forec$.mean[1] + (residuals_std * 1.96)
# Display calculated confidence interval
sprintf(
"Calculated confidence interval: [%f, %f]",
lower_bound,
upper_bound
)## [1] "Calculated confidence interval: [76871.012478, 113502.102384]"The manually calculated interval is nearly identical to the interval produced by R, confirming that both methods provide similar results. The small differences are likely due to rounding or slight variations in how the standard deviation and confidence intervals were calculated.
Question 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.
greece_exports <- global_economy %>% filter(Country == "Greece")
greece_exports %>% autoplot(Exports) +
labs(title = "Exports of Greece", y = "Exports (in USD)", x = "Year")
The data shows an overall upward trend in Greece’s exports from the 1960s to the present, indicating steady growth over time, particularly after 2000. There are noticeable declines in export levels during certain periods, notably in the late 1970s and early 1990s. The drop in the 1970s can be attributed to Greece’s dictatorship from 1967 to 1974, which disrupted the economy. A sharp rise in exports occurs after 2000, with continuous growth throughout the 2010s. This boost is largely due to Greece joining the Eurozone and therefore adopting the euro in 2001, which boosted trade and economic integration with other member countries.
b) Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
fit_ANN <- greece_exports %>%
model(ETS(Exports ~ error("A") + trend("N") + season("N")))
# Generate and plot forecasts
fc_ANN <- forecast(fit_ANN, h = "5 years")
autoplot(fc_ANN) +
labs(title = "Forecasts of Greece Exports (ETS(A,N,N))", y = "Exports (in USD)", x = "Year")
c)Compute the RMSE values for the training data.
rmse_ANN <- accuracy(fit_ANN)$RMSE
cat("The RMSE for the ETS(A,N,N) model is:", rmse_ANN, "\n")## The RMSE for the ETS(A,N,N) model is: 1.816383d) 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_AAN <- greece_exports %>%
model(ETS(Exports ~ error("A") + trend("A") + season("N")))
# Generate and plot forecasts
fc_AAN <- forecast(fit_AAN, h = "5 years")
autoplot(fc_AAN) +
labs(title = "Forecasts of Greece Exports (ETS(A,A,N))", y = "Exports (in USD)", x = "Year")
rmse_AAN <- accuracy(fit_AAN)$RMSE
cat("The RMSE for the ETS(A,A,N) model is:", rmse_AAN, "\n")## The RMSE for the ETS(A,A,N) model is: 1.765767ETS(A,A,N)
Pros: This model’s ability to account for trends can provide more accurate forecasts in datasets where trends are significant. Cons: The added complexity can lead to overfitting, especially if the trend does not persist into the future.
ETS(A,N,N)
Pros: This simpler model is easier to interpret and may be more robust to fluctuations in the data. It might perform adequately when the data lacks a strong trend or when simplicity is preferred for interpretability.Cons: Its inability to capture trends may result in poorer forecasts in datasets where trends are indeed present.
Given the lower RMSE for the ETS(A,A,N) model (1.765767), it is likely the better choice for this dataset, as it effectively captures the upward trend in Greece’s exports. However, if interpretability and simplicity are prioritized, the ETS(A,N,N) model may still be a valid option, particularly if the forecasts do not require that much of precision.
e) Compare the forecasts from both methods. Which do you think is best?
The ETS(A,N,N) forecast shows a flat line, indicating that the model does not account for any trend in the data. The projections remain relatively constant over the forecast horizon, suggesting limited growth in exports. On the other hand the ETS(A,A,N) model shows a clear upward trajectory in the forecast, which aligns better with the long-term trend observed in the data. The forecast suggests a gradual increase in Greece’s exports over time. In conclusion, the ETS(A,A,N) model is preferred for forecasting Greece’s exports due to its ability to incorporate the trend present in the data, leading to more realistic forecasts.
# Generate forecasts for the next periods
forecast_ANN <- forecast(fit_ANN, h = 4) # Forecast for the next 4 periods
forecast_AAN <- forecast(fit_AAN, h = 4)
## Calculated by R
forecast_ANN %>% hilo(95) %>% pull('95%') %>% head(1)## <hilo[1]>
## [1] [29.5966, 36.84272]95forecast_AAN %>% hilo(95) %>% pull('95%') %>% head(1)## <hilo[1]>
## [1] [30.07503, 37.24849]95## Calculated by hand
forecast_AAN <- fit_AAN %>% forecast(h = 4)
forecast_ANN <- fit_ANN %>% forecast(h = 4)
y_hat_AAN <- forecast_AAN$.mean[1]
rmse_AAN <- accuracy(fit_AAN)$RMSE
# Get the first forecast value and RMSE for ETS(A,N,N)
y_hat_ANN <- forecast_ANN$.mean[1]
rmse_ANN <- accuracy(fit_ANN)$RMSE
z <- 1.96
# Calculate prediction intervals for ETS(A,N,N)
lower_ANN <- y_hat_ANN - z * rmse_ANN
upper_ANN <- y_hat_ANN + z * rmse_ANN
# Calculate prediction intervals for ETS(A,A,N)
lower_AAN <- y_hat_AAN - z * rmse_AAN
upper_AAN <- y_hat_AAN + z * rmse_AAN
cat("Prediction Interval for ETS(A,N,N): [", lower_ANN, ", ", upper_ANN, "]\n")## Prediction Interval for ETS(A,N,N): [ 29.65955 , 36.77977 ]cat("Prediction Interval for ETS(A,A,N): [", lower_AAN, ", ", upper_AAN, "]\n")## Prediction Interval for ETS(A,A,N): [ 30.20086 , 37.12266 ]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.
For the ETS(A,N,N) model the manual calculation slightly overestimates the lower bound and underestimates the upper bound compared to the values produced by R. For the ETS(A,N,N) model the manual calculation also shows a higher lower bound and a lower upper bound compared to R. Both methods provide a reasonably close range of estimates, indicating that our manual calculations are on track, but they might not perfectly match the intervals produced by R due to slight differences in the RMSE calculations or the method of interval estimation.
Question 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.
Lambda_china <- global_economy %>%
filter(Country == "China") %>%
features(GDP, features = guerrero) %>%
pull(lambda_guerrero)
China_fit <- global_economy %>%
filter(Country == "China") %>%
model(`Standard` = ETS(GDP ~ error("A") + trend("N") + season("N")),
`Holt's method` = ETS(GDP ~ error("A") + trend("A") + season("N")),
`Damped Holt's method` = ETS(GDP ~ error("A") + trend("Ad", phi = 0.8) + season("N")),
`Box-Cox` = ETS(box_cox(GDP,Lambda_china) ~ error("A") + trend("Ad") + season("N")),
`Damped Box-Cox` = ETS(box_cox(GDP,Lambda_china) ~ error("A") + trend("Ad", phi = 0.8) + season("N")))
China_fc <- China_fit %>%
forecast(h = 20)
China_fc %>%
autoplot(global_economy, level = NULL) +
labs(title="Chinese GDP") +
guides(colour = guide_legend(title = "Forecast"))
In our plot, methods like Damped Box-Cox and Damped Holt’s method provide more conservative GDP growth forecasts than models without damping. The use of a Box-Cox transformation (seen as the red line) leads to higher forecasts, which could imply that the transformation shows faster growth under certain conditions.The use of Holt’s method (shown in blue) typically results in higher growth projections, capturing the underlying linear trend. The Standard ETS method represents the baseline model without any trend dampening or Box-Cox transformation, showing a medium growth.
Question 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?
gas_data <- aus_production %>%
filter(Gas > 0)
aus_production %>%
autoplot(Gas)
# Fit different ETS models including multiplicative seasonality and damped trend
fitted_models <- gas_data %>%
model(
simple_additive = ETS(Gas ~ error('A') + trend("A") + season("N")),
simple_mult = ETS(Gas ~ error("M") + trend("A") + season("N")),
seasonality_add = ETS(Gas ~ error('A') + trend("A") + season("A")),
seasonality_mult = ETS(Gas ~ error("M") + trend("A") + season("M")),
mult_season_damp = ETS(Gas ~ error("M") + trend("Ad") + season("M"))
)
# Generate forecasts for the next 10 time periods
forecasts <- fitted_models %>% forecast(h = 10)
# Plot the forecasts with the original data
forecasts %>%
autoplot(gas_data, level = NULL) +
labs(y = "Gas Production", title = "Forecasts of Australian Gas Production") +
facet_wrap(~ .model, scales = "free_y")
accuracy(fitted_models)## # A tibble: 5 × 10
## .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 simple_additive Traini… 0.442 15.7 11.7 1.05 14.3 2.11 2.07 0.118
## 2 simple_mult Traini… -0.0787 16.4 11.9 -1.32 13.5 2.14 2.16 0.0888
## 3 seasonality_add Traini… 0.00525 4.76 3.35 -4.69 10.9 0.600 0.628 0.0772
## 4 seasonality_mult Traini… -0.115 4.60 3.02 0.199 4.08 0.542 0.606 -0.0131
## 5 mult_season_damp Traini… -0.00439 4.59 3.03 0.326 4.10 0.544 0.606 -0.0217In the seasonality_mult and mult_season_damp models, the multiplicative seasonality proves crucial. This is proved by lower RMSE, MAE, and MAPE values compared to additive models (seasonality_add and simple_additive), indicating that the multiplicative seasonality better captures the seasonal pattern of gas production. The gas production data shows that seasonal changes get bigger as production increases. This pattern is better captured using a model with multiplicative seasonality. The mult_season_damp model has slightly lower RMSE (4.591840) than the seasonality_mult model (4.595113). This suggests that while the trend in gas production continues to grow, it might begin to stabilize in the future, making a damped trend more appropriate for longer-term forecasts. Overall, based on the performance metrics (especially RMSE and MAPE), the mult_season_damp model is the best choice as it combines the strengths of multiplicative seasonality and a damped trend, providing a more accurate and realistic forecast for the gas production data.
Question 8
Recall your retail time series data (from Exercise 7 in Section 2.10).
set.seed(2601)
myseries <- aus_retail %>%
filter(`Series ID` == sample(aus_retail$`Series ID`,1))a) Why is multiplicative seasonality necessary for this series?
myseries %>% autoplot()## Plot variable not specified, automatically selected `.vars = Turnover`
Multiplicative seasonality is necessary here because the seasonal variations grow larger as the overall level of the series increases. This means that the seasonal patterns are proportional to the level of turnover, making a multiplicative model a better fit to capture these changing seasonal effects over time.
b)Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
# Fit the Holt-Winters’ multiplicative models
fit <- myseries %>% model(
Multi_Season = ETS(Turnover ~ error("M") + trend("A") + season("M")),
Multi_Season_Damp = ETS(Turnover ~ error("M") + trend("Ad") + season("M"))
)
# Forecast the next 32 months
fc <- fit %>% forecast(h = 32)
# Plot the forecasts for both models along with the original data
fc %>%
autoplot(myseries, level = NULL) +
labs(y = "AUD") +
theme_minimal()
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 Multi_Season 3.05
## 2 Multi_Season_Damp 3.07The Multiplicative Seasonal Model has a lower RMSE (3.045007) compared to the Multiplicative Seasonal with Damped Trend Model (3.074567). This indicates that the Multiplicative Seasonal Model is a better fit for the data, as it has smaller forecast errors.
d)Check that the residuals from the best method look like white noise.
best_fit <- myseries %>% model(
Multi_Season = ETS(Turnover ~ error("M") + trend("A") + season("M"))
)
gg_tsresiduals(best_fit)
The residuals appear to fluctuate around zero with no clear trend or pattern.In the ACF plot, the bars mostly fall within the confidence bands.The histogram appears roughly centered around zero with no significant skewness, which is another characteristic of normally distributed residuals.
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?
set.seed(2601)
train <- myseries %>%
filter(year(Month) < 2011)
# Test data: from 2011 onwards
test <- myseries %>%
filter(year(Month) >= 2011)
fit <- train %>%
model(
'Holt Winters Multiplicative Method' = ETS(Turnover ~ error('M') + trend('A') + season('M')),
'Holt Winters Damped Method' = ETS(Turnover ~ error('M') + trend('Ad') + season('M')),
'Seasonal Naive' = SNAIVE(Turnover)
)
comp <- anti_join(myseries, train, by = c('State', 'Industry', 'Series ID', 'Month', 'Turnover'))
fc <- fit %>% forecast(comp)
autoplot(comp, Turnover) +
autolayer(fc, level = NULL) +
labs(title = 'Forecast Method Comparison')
# Extract forecasted values for each model and convert to a tibble
forecast_values <- fc %>%
as_tibble() %>%
select(Month, .model, .mean)
# Convert test data to tibble
test_tibble <- test %>%
select(Month, Turnover) %>%
as_tibble()
# Merge forecasted values with actual test values by 'Month'
comparison <- test_tibble %>%
left_join(forecast_values, by = "Month")
# Calculate squared errors between actual and forecasted values
comparison <- comparison %>%
mutate(squared_error = (Turnover - .mean)^2)
# Calculate RMSE manually for each model
rmse_manual <- comparison %>%
group_by(.model) %>%
summarise(RMSE = sqrt(mean(squared_error, na.rm = TRUE)))
rmse_manual## # A tibble: 3 × 2
## .model RMSE
## <chr> <dbl>
## 1 Holt Winters Damped Method 10.2
## 2 Holt Winters Multiplicative Method 13.1
## 3 Seasonal Naive 6.12The Seasonal Naive model has the lowest RMSE at 6.12. I could not beat the Seasonal Naive approach from Exercise 7 in Section 5.11. The Seasonal Naive model outperformed the Holt Winters models in terms of RMSE.
Question 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_train <- train %>%
features(Turnover, features = guerrero) %>%
pull(lambda_guerrero)
box_cox_myseries <- train %>%
mutate(
bc = box_cox(Turnover, lambda_train)
)
fit2 <- box_cox_myseries %>%
model(
'STL Box-Cox' = STL(bc ~ season(window = 'periodic'), robust = TRUE),
'ETS Box-Cox' = ETS(bc)
)
fit3 <- box_cox_myseries %>%
model(
'Holt Winters Multiplicative Method' = ETS(Turnover ~ error('M') + trend('A') + season('M'))
)
accuracy(fit2)## # A tibble: 2 × 12
## State Industry .model .type ME RMSE MAE MPE MAPE MASE RMSSE
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Wester… Footwea… STL B… Trai… 0.00329 0.0562 0.0426 0.0570 1.24 0.354 0.375
## 2 Wester… Footwea… ETS B… Trai… -0.00123 0.0676 0.0529 -0.0545 1.52 0.440 0.451
## # ℹ 1 more variable: ACF1 <dbl>accuracy(fit3)## # A tibble: 1 × 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 West… Footwea… Holt … Trai… 0.175 2.68 1.72 0.00277 4.63 0.473 0.500 0.0586