Data624_HW05
Mahmud Hasan Al Raji
2024-10-06
In this assignment, the problems 8.1, 8.5, 8.6, 8.7, 8.8, and 8.9 have been solved from the Hyndman book. (book link:https://otexts.com/fpp3/toolbox-exercises.html)
Libraries
library(fpp3)## Registered S3 method overwritten by 'tsibble':
## method from
## as_tibble.grouped_df dplyr## ── Attaching packages ──────────────────────────────────────────── fpp3 1.0.0 ──## ✔ tibble 3.2.1 ✔ tsibble 1.1.5
## ✔ dplyr 1.1.4 ✔ tsibbledata 0.4.1
## ✔ tidyr 1.3.1 ✔ feasts 0.3.2
## ✔ lubridate 1.9.3 ✔ fable 0.3.4
## ✔ ggplot2 3.5.1 ✔ fabletools 0.4.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()Question-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 ℓ0,and generate forecasts for the next four months.
#head(aus_livestock)
aus_pigs_slaughter<-aus_livestock|>filter(Animal=="Pigs",State=="Victoria")
#aus_pigs_slaughter
# Estimate parameters
fit <- aus_pigs_slaughter |>
model(ETS(Count ~ error("A") + trend("N") + season("N")))
# Get model fit report
report(fit)## 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.07# Generate forecast
fc <- fit |>
forecast(h = 4)The optimal parameter estimates α =0.32 and ℓ0= 1.01e+05.
# Plot forecasted values along with actual observed values
fc |>
autoplot(aus_pigs_slaughter) +
geom_line(aes(y = .fitted), col="#D55E00",
data = augment(fit)) +
labs(y="Slaughter Count", title="The number of pigs slaughtered in Victoria") +
guides(colour = "none")
#Plot forecasted values from 2015 January onward along with actual observed values
fc|> autoplot(aus_pigs_slaughter |> filter(Month >= yearmonth("2015 Jan"))) +
geom_line(aes(y = .fitted, col="#D55E00"),
data = augment(fit) |> filter(Month >= yearmonth("2015 Jan"))) +
labs(y = "Slaughter Count",
title = "The number of pigs slaughtered in Victoria (from Jan 2015)") +
guides(colour = "none")
The forecasts for the monthly pig slaughter counts in Victory are plotted for the period from January 2019 to April 2019. The orange colored line represents the one-step-ahead fitted values in both the plots. The second plot is plotted with same fitted values for smaller time period just to see the forecasts more clearly.The colored/shaded area around the forecasts reflects an 80% probability that almost all the actual values will fall within the prediction intervals of the forecasts.
- 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.
Manual computation of 95% prediction intervals for the first forecast:
#Forecast for 4 months
fc <- fit |>
forecast(h =4)
#fc
#Get mean
y_hat <- fc$.mean[1]
#Get residuals
#augment(fit)
augment_fit <- augment(fit)
#Get standard dev using residuals
s <- sd(augment_fit$.resid)
#Calculate 95% prediction intervals
upper_95_level <- y_hat + (s * 1.96)
lower_95_level <- y_hat - (s * 1.96)
#Lower 95 interval
lower_95_level## [1] 76871.01#Upper interval
upper_95_level## [1] 113502.195% confidence intervals for the first first forecast produced by R:
# Extract lower and upper 95% prediction intervals for first forecast
first_forecast_95<-fc %>% hilo() %>% pull(`95%`) %>% .[1]
first_forecast_95## <hilo[1]>
## [1] [76854.79, 113518.3]95The manually computed 95% prediction interval for the first forecast is [76871.01, 113502.1], while the interval generated by R is [76854.79, 113518.3]. Both intervals are nearly identical with only minimal differences. This shows that the manual calculation is consistent with R generated values.
Question-8.5:
Data set global_economy contains the annual Exports from many countries. Select one country to analyse.
#?global_economy
bd_economy <- global_economy|>
filter(Country == 'Bangladesh')
#sum(is.na(bd_economy))- Plot the Exports series and discuss the main features of the data.
#head(bd_economy)
bd_economy|> autoplot(Exports)+
labs(y = "Exports (% of GDP)", title = "Exports of goods and services of Bangladesh")
The time series above represents Bangladesh’s exports of goods and services as a percentage of its GDP from 1960 to 2017. It shows a decline in exports from 1960 to 1975, followed by a rise in exports until 2012, with four noticeable dips during this period. After 2012, exports consistently decline through 2017.
- Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
#Estimate parameters
fit_bd <- bd_economy|>
model(ETS(Exports ~ error("A") + trend("N") + season("N")))
#Create forecast
fc <- fit_bd |>
forecast(h = 4)
#Plot forecasts
fc |>
autoplot(bd_economy) +
labs(y="Exports (% of GDP)", title="Exports of goods and services of Bangladesh") +
guides(colour = "none")
- Compute the RMSE values for the training data.
fit_bd |> 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 Bangladesh "ETS(Expo… Trai… 0.0869 1.25 0.945 -0.891 12.0 0.983 0.991 0.0148Train_RMSE <- round((fit_bd |> accuracy())$RMSE,2)
print(paste("The RMSE value of the ETS(A,N,N) model is ", Train_RMSE))## [1] "The RMSE value of the ETS(A,N,N) model is 1.25"- 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.
comparison_fit <- bd_economy |>
model(
ANN = ETS (Exports ~ error("A") + trend("N") + season("N")),
AAN = ETS (Exports ~ error("A") + trend("A") + season("N"))
)
accuracy(comparison_fit)## # 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 Bangladesh ANN Training 0.0869 1.25 0.945 -0.891 12.0 0.983 0.991 0.0148
## 2 Bangladesh AAN Training 0.00557 1.25 0.954 -1.90 12.2 0.993 0.989 0.0135The AAN model demonstrates a better performing model compared to the ANN model in terms of RMSE. A lower RMSE indicates that the AAN model’s predictions are closer to the actual values. It basically reflects that the AAN model has a more accurate forecasting capability.
In terms of merits, in general, the ANN model is straightforward and easy to implement. It does not capture complex patterns effectively. On the other hand, the AAN model has increased flexibility and adaptability. It can capture non-linear trends and seasonal variations. In the context of this dataset, there is no pattern in the time series observed, though. But the time series shows a non-linear upward trend with irregular ups and downs. So, here I think AAN will be a good choice to predict future values. The lower RMSE value also confirms that the AAN model is likely the more robust choice here for accurate forecasting.
- Compare the forecasts from both methods. Which do you think is best?
comparison_fit |>
forecast(h=4) |>
autoplot(bd_economy, level=NULL) +
labs(title="Forecasts Comparison",
subtitle = "Exports of Bangladesh")
Comparing the nature of the time series it seems to me that the forecasts from the AAN model will perform better than those produced by the ANN model. Because, the time series here shows a non-linear upward trend with irregular ups and downs.The forecasts from the AAN model also show a upward trend, which seems rational with the past pattern of downward trends followed by upward movement, reflecting the overall upward direction of the series. On the other hand, the forecasts from the ANN show a flat line, which does not seem to properly reflect the nature of this time series.
- 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.
Manual computation of 95% prediction intervals of first forecast produced by ANN model:
# Fit ANN model
fit_bd_ANN <- bd_economy|>
model(ETS(Exports ~ error("A") + trend("N") + season("N")))
#Forecast for 4 months
fc_ANN <- fit_bd_ANN |>
forecast(h =4)
#fc
#Get mean
y_hat_ANN <- fc_ANN$.mean[1]
#Get residuals
#augment(fit)
augment_fit_ANN <- augment(fit_bd_ANN)
#Get standard dev using residuals
#s <- sd(augment_fit_ANN$.resid)
#Get RMSE from residuals
rmse1 <- sqrt(mean(augment_fit_ANN$.resid^2))
#Calculate 95% prediction intervals
upper_95_level_ANN <- y_hat_ANN + (rmse1 * 1.96)
lower_95_level_ANN <- y_hat_ANN - (rmse1 * 1.96)
#Lower 95 interval
lower_95_level_ANN## [1] 12.58008#Upper interval
upper_95_level_ANN## [1] 17.49246Calculate 95% confidence intervals for the first first forecast produced by R for ANN model:
# Extract lower and upper 95% prediction intervals for first forecast
first_forecast_95_ANN<-fc_ANN %>% hilo() %>% pull(`95%`) %>% .[1]
first_forecast_95_ANN## <hilo[1]>
## [1] [12.53665, 17.53589]95The manually computed 95% prediction intervals for the first forecast produced by the ANN model are [12.52301, 17.54953], while the R-generated intervals are [12.44063, 17.63191]. Both intervals are very close, with only slight differences in the bounds.The R-generated interval is slightly wider. But overall, both approaches produce similar results, indicating consistency between the manual calculation and R generated values.
Manual computation of 95% prediction intervals of first forecast produced by AAN model:
# Fit ANN model
fit_bd_AAN <- bd_economy|>
model(ETS(Exports ~ error("A") + trend("A") + season("N")))
#Forecast for 4 months
fc_AAN <- fit_bd_AAN |>
forecast(h =4)
#fc
#Get mean
y_hat_AAN <- fc_AAN$.mean[1]
#Get residuals
#augment(fit)
augment_fit_AAN <- augment(fit_bd_AAN)
#Get standard dev using residuals
#s <- sd(augment_fit_AAN$.resid)
#Get RMSE from residuals
rmse2 <- sqrt(mean(augment_fit_AAN$.resid^2))
#Calculate 95% prediction intervals
upper_95_level_AAN <- y_hat_AAN + (rmse2 * 1.96)
lower_95_level_AAN <- y_hat_AAN - (rmse2 * 1.96)
#Lower 95 interval
lower_95_level_AAN## [1] 12.6639#Upper interval
upper_95_level_AAN## [1] 17.56622Calculate 95% confidence intervals for the first first forecast produced by R for ANN model:
# Extract lower and upper 95% prediction intervals for first forecast
first_forecast_95_AAN<-fc_AAN %>% hilo() %>% pull(`95%`) %>% .[1]
first_forecast_95_AAN## <hilo[1]>
## [1] [12.57479, 17.65533]95The manually computed 95% prediction intervals for the first forecast produced by the AAN model are [12.86225, 17.77228], while the R-generated intervals are [12.69279, 17.94174]. Both intervals are very close, with only slight differences in the bounds.Again, the R-generated interval is slightly wider. Again, both approaches produce similar results, indicating consistency between the manual calculation and R generated values.
Question-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.]
Generate time series of China’s GDP:
# head(global_economy)
china_gdp <- global_economy|>
filter(Country == "China")
china_gdp |> autoplot(GDP) +
labs(title="GDP of China")
The time series of China’s GDP shows an upward trend. The time series does not reflect any seasonality in the data. I will try to understand the differences of forecasts produced by various models utilizing ETS() functions to see how much the forecasts change with damped trend, or with a Box-Cox transformation.
#Get lambda
lambda <- global_economy|> filter(Country == "China") |> features(GDP, features = guerrero)|> pull(lambda_guerrero)
#Fit different ETS models
fit_china_gdp <- china_gdp|>
model(`Simple` = 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) ~ error("A") + trend("A") + season("N")),
`Box-Cox Damped` = ETS(box_cox(GDP,lambda) ~ error("A") + trend("Ad", phi = 0.8) + season("N")))
#Produce models' forecasts
fc_china_gdp <- fit_china_gdp|>forecast(h = 10)
# Plot forecasts
fc_china_gdp |> autoplot(china_gdp, level = NULL) +
labs(title="China's GDP Forecasts") + guides(colour = guide_legend(title = "GDP Forecast")) 
# Plot forecasts faceted by models
fc_china_gdp |> autoplot(china_gdp, level = NULL) +
labs(title="China's GDP Forecasts") + guides(colour = guide_legend(title = "GDP Forecast"))+facet_wrap(~.model, scales = "free_y") 
Observing the nature of the time series curve, it seems that the simple exponential smoothing model is failing to produce accurate forecasts for this dataset. The Box-Cox method is generating overly high forecasts. The Damped Holt’s method is producing slightly curved forecasts that miss the past upward trend a bit. The forecasts produced by Box-Cox Damped and Holt’s method are appearing almost similar, with the Box-Cox Damped line showing a slightly steeper incline. I would say that forecasts from either Box-Cox Damped or Holt’s method are likely to perform well in this case.However, analyzing the model’s performance accuracy will give us the more appropriate result.
Question-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?
Plot time series of quarterly Autralian gas production:
#head(aus_production)
#?aus_production
aus_production |> autoplot(Gas)+
labs(title=" Quarterly Gas Production of Australia")
The most noticeable aspect of this data is the change in the seasonal trend. As time goes on, the seasonal variation with increasing trend becomes more evident. Because of this, a multiplicative seasonal component is necessary for the ETS model in this case.
Apply Holt-Winters’ method with multiplicative seasonality, and with a damped trend and multiplicative seasonality:
fit_gas <- aus_production |>
model(
Multiplicative = ETS(Gas ~ error("M") + trend("A") + season("M")),
`Damped Multiplicative` = ETS(Gas ~ error("M") + trend("Ad") + season("M"))
)
fc_gas <- fit_gas |> forecast(h = "5 years")
fc_gas |>
autoplot(aus_production, level = NULL) +
labs(title=" Quarterly Gas Production of Australia") +
guides(colour = guide_legend(title = "Forecast"))
Get models’ performance metrics:
report(fit_gas)## Warning in report.mdl_df(fit_gas): Model reporting is only supported for
## individual models, so a glance will be shown. To see the report for a specific
## model, use `select()` and `filter()` to identify a single model.## # A tibble: 2 × 9
## .model sigma2 log_lik AIC AICc BIC MSE AMSE MAE
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Multiplicative 0.00324 -831. 1681. 1682. 1711. 21.1 32.2 0.0413
## 2 Damped Multiplicative 0.00329 -832. 1684. 1685. 1718. 21.1 32.0 0.0417Based on the AIC values, it seems that the Multiplicative Holt-Winters’ method performs slightly better than the Holt-Winters’ damped trend method.Therefore, the simpler standard Multiplicative Holt-Winters’ method is recommended here.
Question-8.8
Recall your retail time series data (from Exercise 7 in Section 2.10).
set.seed(12345)
myseries <- aus_retail |>
filter(`Series ID` == sample(aus_retail$`Series ID`,1))
#head(myseries)
myseries |> autoplot(Turnover)+
labs(title=" Monthly Turnover of Australia's Retail Market")
- Why is multiplicative seasonality necessary for this series?
The time series for the Australian retail market possesses changes in the seasonal trend. As time goes on, the seasonal variation with increasing trend becomes more evident. Because of this, a multiplicative seasonal component is necessary for this series.
- Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
fit_myseries <- myseries |>
model(Multiplicative = ETS(Turnover ~ error("M") + trend("A") + season("M")),
`Damped Multiplicative` = ETS(Turnover ~ error("M") + trend("Ad") + season("M")))
fit_myseries |>
forecast(h=48) |>
autoplot(myseries, level = NULL) +
labs(title="Monthly Turnover of Australia's Retail Market") +
guides(colour = guide_legend(title = "Forecast"))
The Holt-Winters’ damped trend method seems to be more suitable for producing forecasts for this data, given the seasonal patterns in the series. However, it is difficult here to distinguish between the forecasts produced by the damped and un-damped methods,though.
- Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
accuracy(fit_myseries) |> dplyr::select(.model, RMSE)## # A tibble: 2 × 2
## .model RMSE
## <chr> <dbl>
## 1 Multiplicative 5.59
## 2 Damped Multiplicative 5.63The RMSE value is lower for the simpler Holt-Winters’ multiplicative method. Hence,the un-damped multiplicative is the preferred method here for the one-step forecasts.
- Check that the residuals from the best method look like white noise.
myseries |> model(Multiplicative = ETS(Turnover ~ error("M") + trend("A") + season("M")))|>
gg_tsresiduals() + ggtitle(" Residuals Plots of Turnover Using Multiplicative Method")
The spread of the residuals does not seem constant around the mean. The distribution of residuals seems very close to normal. The autocorrelation function (ACF) plot shows that more than 5% (5/26*100 ~ 19.23%) of the spikes cross the significance bounds, it suggests that residuals are not uncorrelated. Considering these matters, it can be said the residuals don’t look like white noise.
- 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?
#Training data
myseries_train <- myseries |>
filter(year(Month) < 2011)
#Test data
myseries_test <- myseries |>
filter(year(Month) >= 2011)
#Fit SNAIVE and Multiplicative models
fit_myseries_train <- myseries_train |>
model(Multiplicative = ETS(Turnover ~ error("M") + trend("A") + season("M")),
SNAIVE = SNAIVE(Turnover))
#Produce forecasts
fc_myseries_train <- fit_myseries_train |>
forecast(h=15)
fc_myseries_train |> autoplot(myseries_train, level = NULL)+labs(title="Turnover Forecasts")+ guides(colour = guide_legend(title = "Forecasts"))
Get model’s performance metric RMSE for training data:
accuracy(fit_myseries_train) |>
dplyr::select(.type, .model, RMSE)## # A tibble: 2 × 3
## .type .model RMSE
## <chr> <chr> <dbl>
## 1 Training Multiplicative 5.16
## 2 Training SNAIVE 11.9Get model’s performance metric RMSE on test data:
fc_myseries_train |> accuracy(myseries_test) |>
dplyr::select(.type, .model, RMSE)## # A tibble: 2 × 3
## .type .model RMSE
## <chr> <chr> <dbl>
## 1 Test Multiplicative 15.6
## 2 Test SNAIVE 10.9Based on the RMSE values, the Multiplicative model performs well during training but does not outperform the SNAIVE method on the test data.
Question-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?
#Get lambda
lambda1<- myseries_train |>
features(Turnover, features = guerrero) |>
pull(lambda_guerrero)
#Box-cox transformation to the training data set
myseries_train_boxcox <- myseries_train
myseries_train_boxcox$Turnover <- box_cox(myseries_train$Turnover,lambda1)
#Box-cox transformation to test data set
myseries_test_boxcox<-myseries_test
myseries_test_boxcox$Turnover <- box_cox(myseries_test$Turnover,lambda1)
#Apply STL decomposition
stl_decomp <- myseries_train_boxcox |>
model(stl = STL(Turnover))
components(stl_decomp) %>% autoplot()
The plot above shows that applying the Box-Cox transformation to the series has resulted in a nearly constant seasonal effect, as well as a reduction in variability in the remainder.
#Fit multiplicative model
fit_train_boxcox <- myseries_train_boxcox |>
model(
"Multiplicative Method" = ETS(Turnover ~ error('M') + trend('A') + season('M'))
)
fc_train_boxcox <- fit_train_boxcox |>
forecast(h = 15)
fc_train_boxcox |> autoplot(myseries_train_boxcox, level = NULL) +
labs(title = "Turnover Forecasts") +
guides(colour = guide_legend(title = "Forecasts"))
See model performance metric ‘RMSE’ on model training :
accuracy(fit_train_boxcox) |> dplyr::select(.type, .model, RMSE)## # A tibble: 1 × 3
## .type .model RMSE
## <chr> <chr> <dbl>
## 1 Training Multiplicative Method 0.0805Applying the Box-Cox transformation to the time series data has enhanced the multiplicative model’s performance by significantly reducing the RMSE value. The RMSE previously recorded was 5.16 for the model with the training dataset, but it has now decreased to 0.081.
See model performance metric ‘RMSE’ on test data :
fc_train_boxcox |> accuracy(myseries_test_boxcox) |>
dplyr::select(.type, .model, RMSE)## # A tibble: 1 × 3
## .type .model RMSE
## <chr> <chr> <dbl>
## 1 Test Multiplicative Method 0.197Previously, the SNAIVE method demonstrated better performance on the test data set, achieving an RMSE of 10.863. In contrast, the Multiplicative model had a higher RMSE of 15.633. However, after applying the Box-Cox transformation, the Multiplicative model has become the best performer on the test data set. It has achieved a significantly lower RMSE of 0.197.