DATA 624 HW 5
Rebecca Bronstein
2026-03-08
8.8 Exercises
library(fpp3)## Registered S3 method overwritten by 'tsibble':
## method from
## as_tibble.grouped_df dplyr## ── Attaching packages ──────────────────────────────────────────── fpp3 1.0.2 ──## ✔ tibble 3.3.1 ✔ tsibble 1.1.6
## ✔ dplyr 1.1.4 ✔ tsibbledata 0.4.1
## ✔ tidyr 1.3.2 ✔ feasts 0.4.2
## ✔ lubridate 1.9.4 ✔ fable 0.5.0
## ✔ ggplot2 4.0.1## ── 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()Consider the the number of pigs slaughtered in Victoria, available in the
aus_livestockdataset.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.The simple exponential smoothing model was estimated using the ETS A N N specification. The estimated smoothing parameter alpha is 0.3221247 and the estimated initial level l0 is 100646.6. Forecasts for the next four months are constant because simple exponential smoothing assumes no trend or seasonality, so all forecasts equal the final estimated level.
# Extract Victoria pig slaughter data vic_pigs <- aus_livestock %>% filter(Animal == "Pigs", State == "Victoria") # Fit simple exponential smoothing fit <- vic_pigs %>% model(SES = ETS(Count ~ error("A") + trend("N") + season("N"))) 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#next 4 months fc <- fit %>% forecast(h = "4 months") autoplot(fc, vic_pigs)
- 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.
#calculate the 95% prediction interval hilo(fc, 95)## # A tsibble: 4 x 7 [1M] ## # Key: Animal, State, .model [1] ## Animal State .model Month ## <fct> <fct> <chr> <mth> ## 1 Pigs Victoria SES 2019 Jan ## 2 Pigs Victoria SES 2019 Feb ## 3 Pigs Victoria SES 2019 Mar ## 4 Pigs Victoria SES 2019 Apr ## # ℹ 3 more variables: Count <dist>, .mean <dbl>, `95%` <hilo>The residual variance is 87480760 which gives a residual standard deviation s of 9353.6. Using the approximation yhat plus or minus 1.96 times s the 95 percent prediction interval for the first forecast is approximately 82567 to 119233. This interval is very close to the interval produced by R with small differences occurring because the ETS model calculates prediction intervals using the full state space variance rather than the simple approximation.
Write your own function to implement simple exponential smoothing. The function should take arguments
y(the time series),alpha(the smoothing parameter α) andlevel(the initial level ℓ0). It should return the forecast of the next observation in the series. Does it give the same forecast asETS()?#model ses_forecast <- function(y, alpha, level) { l <- level for (i in 1:length(y)) { l <- alpha * y[i] + (1 - alpha) * l } return(l) } #applied model alpha <- 0.3221247 level <- 100646.6 ses_forecast(vic_pigs$Count, alpha, level)## [1] 95186.56#compared fit |> forecast(h = 1)## # A fable: 1 x 6 [1M] ## # Key: Animal, State, .model [1] ## Animal State .model Month ## <fct> <fct> <chr> <mth> ## 1 Pigs Victoria SES 2019 Jan ## # ℹ 2 more variables: Count <dist>, .mean <dbl>The model applies the recursive simple exponential smoothing formula to update the level for each observation in the series. The forecast for the next observation is the final level estimate. When the estimated values of alpha = 0.3221247 and the initial level = 100646.6 are used, the function produces a forecast of about 95186.56. This is the same forecast produced by the ETS function, confirming that the model implementation is correct.
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 theETS()function?#modified ses_sse <- function(par, y) { alpha <- par[1] level <- par[2] l <- level sse <- 0 for (i in 1:length(y)) { yhat <- l error <- y[i] - yhat sse <- sse + error^2 l <- alpha * y[i] + (1 - alpha) * l } return(sse) } #optimal y <- vic_pigs$Count opt <- optim( par = c(0.3, mean(y)), fn = ses_sse, y = y, method = "L-BFGS-B", lower = c(0, 0), upper = c(1, Inf) ) opt$par## [1] 3.219965e-01 1.005328e+05The function was modified to calculate the sum of squared errors by computing the one step ahead forecast error at each time step and adding the squared errors together. The optim function was then used to find the values of alpha and the initial level that minimize the SSE. The estimated parameters were alpha = 0.3219965 and l0 = 100532.8. These values are very close to the estimates obtained from the ETS function, which were alpha = 0.3221247 and l0 = 100646.6. The slight differences occur because optim finds the parameters by numerically minimizing the sum of squared errors, while ETS estimates the parameters using a state space maximum likelihood approach.
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.
This function combines the previous two functions by first using optim to estimate the values of alpha and the initial level that minimize the sum of squared errors. After finding the optimal parameters, the function applies the simple exponential smoothing formula to update the level through the series and compute the forecast for the next observation. The function returns the estimated alpha, the initial level, and the forecast.
#combined ses_model <- function(y) { # Function to compute SSE ses_sse <- function(par, y) { alpha <- par[1] level <- par[2] l <- level sse <- 0 for (i in 1:length(y)) { yhat <- l error <- y[i] - yhat sse <- sse + error^2 l <- alpha * y[i] + (1 - alpha) * l } return(sse) } # Optimize alpha and initial level opt <- optim( par = c(0.3, mean(y)), fn = ses_sse, y = y, method = "L-BFGS-B", lower = c(0, 0), upper = c(1, Inf) ) alpha <- opt$par[1] level <- opt$par[2] # Compute final level for forecast l <- level for (i in 1:length(y)) { l <- alpha * y[i] + (1 - alpha) * l } forecast <- l return(list(alpha = alpha, l0 = level, forecast = forecast)) } ses_model(vic_pigs$Count)## $alpha ## [1] 0.3219965 ## ## $l0 ## [1] 100532.8 ## ## $forecast ## [1] 95186.74Data set
global_economycontains the annual Exports from many countries. Select one country to analyse.r- Plot the Exports series and discuss the main features of the data.
# Select China
china_exports <- global_economy |>
filter(Country == "China")
# Plot exports over time
china_exports |>
autoplot(Exports) +
labs(title = "Exports of China Over Time",
y = "Exports (% of GDP)",
x = "Year")
The export series for China shows a strong upward trend over time, reflecting rapid economic growth and increased participation in global trade. The most significant growth occurs from the 1990s to the early 2000s, with some fluctuations caused by global economic conditions such as the 2008 Global Financial Crisis. Since the data is annual, no seasonal pattern is present, but the series clearly demonstrates long-term structural growth.
- Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
# Filter the data for China
china_exports <- global_economy |>
filter(Country == "China")
# Fit ETS(A,N,N) model
fit <- china_exports |>
model(ETS(Exports ~ error("A") + trend("N") + season("N")))
# Forecast the next 10 years
fc <- fit |>
forecast(h = 10)
report(fit)## Series: Exports
## Model: ETS(A,N,N)
## Smoothing parameters:
## alpha = 0.9998999
##
## Initial states:
## l[0]
## 4.308241
##
## sigma^2: 3.7396
##
## AIC AICc BIC
## 315.9713 316.4157 322.1526# Plot the data and forecasts
autoplot(china_exports, Exports) +
autolayer(fc, level = NULL) +
labs(title = "Forecast of China's Exports using ETS(A,N,N)",
y = "Exports (% of GDP)",
x = "Year")
The ETS ANN model was fitted to the exports data for China. The estimated alpha value is 0.9999, which means the model gives almost all the weight to the most recent observation when updating the level. The forecasts remain roughly constant around the most recent export level of about 20 percent of GDP because the model does not include a trend or seasonal component.
Compute the RMSE values for the training data.
#fit accuracy accuracy(fit)## # 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 China "ETS(Exports ~ … Trai… 0.266 1.90 1.26 1.84 9.34 0.983 0.991 0.288#value accuracy(fit) |> select(RMSE)## # A tibble: 1 × 1 ## RMSE ## <dbl> ## 1 1.90The RMSE for the training data was calculated using the fitted ETS(A,N,N) model for China. The RMSE measures the average size of the prediction errors in the same units as the export data, with lower values indicating a better fit of the model to the observed data.
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 ETS(A,A,N) model
fit_trend <- china_exports |>
model(ETS(Exports ~ error("A") + trend("A") + season("N")))
# Compute RMSE for both models
accuracy(fit)## # 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 China "ETS(Exports ~ … Trai… 0.266 1.90 1.26 1.84 9.34 0.983 0.991 0.288accuracy(fit_trend)## # 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 China "ETS(Exports… Trai… -0.0854 1.91 1.25 -0.169 9.57 0.973 0.995 0.232The ETS ANN model assumes no trend, while the ETS AAN model includes a trend component. Because the exports data for China shows periods of growth, the ETS AAN model may fit the data better and produce a lower RMSE. However, the AAN model uses one more parameter, so if the improvement in RMSE is small, the simpler ANN model may still be preferred.
- Compare the forecasts from both methods. Which do you think is best?
The ETS ANN forecast remains mostly flat because the model does not include a trend. The ETS AAN forecast continues the upward movement in the data because it includes a trend component. Since the exports data for China shows growth over time, the ETS AAN model likely produces more realistic forecasts and is the better method for this data.
- 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.
# ETS ANN model
fit <- china_exports |>
model(ETS(Exports ~ error("A") + trend("N") + season("N")))
# ETS AAN model
fit_trend <- china_exports |>
model(ETS(Exports ~ error("A") + trend("A") + season("N")))
# Forecast first step
fc_ann <- forecast(fit, h = 1)
fc_aan <- forecast(fit_trend, h = 1)
fc_ann## # A fable: 1 x 5 [1Y]
## # Key: Country, .model [1]
## Country .model Year
## <fct> <chr> <dbl>
## 1 China "ETS(Exports ~ error(\"A\") + trend(\"N\") + season(\"N\"))" 2018
## # ℹ 2 more variables: Exports <dist>, .mean <dbl>fc_aan## # A fable: 1 x 5 [1Y]
## # Key: Country, .model [1]
## Country .model Year
## <fct> <chr> <dbl>
## 1 China "ETS(Exports ~ error(\"A\") + trend(\"A\") + season(\"N\"))" 2018
## # ℹ 2 more variables: Exports <dist>, .mean <dbl>The 95 percent prediction intervals calculated using RMSE are approximately 12.5 to 27.0 for the ETS ANN model and 11.7 to 27.0 for the ETS AAN model. These intervals are very similar to those produced by R, with small differences due to rounding and the model’s internal calculations. Both methods produce similar uncertainty ranges for the first forecast.
Forecast the Chinese GDP from the
global_economydata set using an ETS model. Experiment with the various options in theETS()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
hwhen forecasting, so you can clearly see the differences between the various options when plotting the forecasts.]
# Filter China GDP data
china_gdp <- global_economy |>
filter(Country == "China")
# Basic ETS model
fit1 <- china_gdp |>
model(ETS(GDP))
# ETS with damped trend
fit2 <- china_gdp |>
model(ETS(GDP ~ trend("Ad", phi = 0.9)))
# ETS with Box Cox transformation
fit3 <- china_gdp |>
model(ETS(box_cox(GDP, 0.3)))
# Forecast far ahead to see differences
fc1 <- forecast(fit1, h = 20)
fc2 <- forecast(fit2, h = 20)
fc3 <- forecast(fit3, h = 20)
# Plot forecasts
autoplot(china_gdp, GDP) +
autolayer(fc1, series = "Standard ETS") +
autolayer(fc2, series = "Damped trend") +
autolayer(fc3, series = "Box Cox")## Warning in ggdist::geom_lineribbon(without(intvl_mapping, "colour_ramp"), :
## Ignoring unknown parameters: `series`## Warning in geom_line(mapping = without(mapping, "shape"), data =
## unpack_data(object[single_row[["FALSE"]], : Ignoring unknown parameters:
## `series`## Warning in ggdist::geom_lineribbon(without(intvl_mapping, "colour_ramp"), :
## Ignoring unknown parameters: `series`## Warning in geom_line(mapping = without(mapping, "shape"), data =
## unpack_data(object[single_row[["FALSE"]], : Ignoring unknown parameters:
## `series`## Scale for fill_ramp is already present.
## Adding another scale for fill_ramp, which will replace the existing scale.## Warning in ggdist::geom_lineribbon(without(intvl_mapping, "colour_ramp"), : Ignoring unknown parameters: `series`
## Ignoring unknown parameters: `series`## Scale for fill_ramp is already present.
## Adding another scale for fill_ramp, which will replace the existing scale.
Forecasts of GDP for China were produced using several ETS models. The damped trend model slows the growth rate in long term forecasts, while the Box Cox transformation stabilizes the variance and smooths large increases. These options change how quickly the forecasts grow over time and help control unrealistic long term trends.
- Find an ETS model for the Gas data from
aus_productionand forecast the next few years. Why is multiplicative seasonality necessary here? Experiment with making the trend damped. Does it improve the forecasts?
#ETS Fit and Forecast
# Plot the data
aus_production |>
autoplot(Gas)
# Fit ETS model
gas_fit <- aus_production |>
model(ETS(Gas))
report(gas_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.389# Forecast next 10 years (40 quarters)
gas_fc <- gas_fit |>
forecast(h = 40)
# Plot forecasts
gas_fc |>
autoplot(aus_production)
#Dampedge Trend
gas_damped <- aus_production |>
model(ETS(Gas ~ error("M") + trend("Ad") + season("M")))
report(gas_damped)## Series: Gas
## Model: ETS(M,Ad,M)
## Smoothing parameters:
## alpha = 0.6489044
## beta = 0.1551275
## gamma = 0.09369372
## phi = 0.98
##
## Initial states:
## l[0] b[0] s[0] s[-1] s[-2] s[-3]
## 5.858941 0.09944006 0.9281912 1.177903 1.07678 0.8171255
##
## sigma^2: 0.0033
##
## AIC AICc BIC
## 1684.028 1685.091 1717.873gas_damped_fc <- gas_damped |>
forecast(h = 40)
gas_damped_fc |>
autoplot(aus_production)
An ETS model was fitted to the Gas series from the aus_production dataset and the selected model was ETS M Ad M. This model includes multiplicative errors, a damped additive trend, and multiplicative seasonality. Multiplicative seasonality is necessary because the seasonal fluctuations increase as the level of the series increases, which can be seen in the plot where peaks and troughs become larger over time. The damped trend slows the long term growth of the forecasts rather than allowing the trend to increase indefinitely. The model with the damped trend also has a slightly lower AIC than the non damped version, suggesting a small improvement in model fit.
Recall your retail time series data (from Exercise 7 in Section 2.10).
Why is multiplicative seasonality necessary for this series?
Multiplicative seasonality is needed because the seasonal fluctuations increase as the overall level of the series grows. Early in the series, the seasonal peaks and drops are smaller, but they become larger as retail sales rise over time. This shows that the seasonal pattern changes in proportion to the level of the data rather than staying constant. A multiplicative seasonal model accounts for this by allowing the seasonal effect to scale with the series level.
Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
library(fpp3) set.seed(12345678) myseries <- aus_retail |> filter(`Series ID` == sample(aus_retail$`Series ID`, 1)) # Plot the series myseries |> autoplot(Turnover)
# Fit additive and multiplicative seasonal models fit <- myseries |> model( additive = ETS(Turnover ~ error("A") + trend("A") + season("A")), multiplicative = ETS(Turnover ~ error("M") + trend("A") + season("M")) ) # Compare model accuracy accuracy(fit)## # 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 Northern T… Clothin… addit… Trai… -0.0223 0.676 0.514 -0.515 6.65 0.586 0.584 ## 2 Northern T… Clothin… multi… Trai… -0.0128 0.613 0.450 -0.469 5.15 0.513 0.529 ## # ℹ 1 more variable: ACF1 <dbl>The above plot shows that seasonal fluctuations grow as the level of retail turnover increases. Because the size of the seasonal pattern scales with the level of the data, a multiplicative seasonal model fits better than an additive one, which assumes constant seasonal variation.
- Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
fit <- myseries |>
model(
hw = ETS(Turnover ~ error("M") + trend("A") + season("M")),
hw_damped = ETS(Turnover ~ error("M") + trend("Ad") + season("M"))
)
accuracy(fit)## # 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 Northern … Clothin… hw Trai… -0.0128 0.613 0.450 -0.469 5.15 0.513 0.529
## 2 Northern … Clothin… hw_da… Trai… 0.0398 0.616 0.444 -0.0723 5.06 0.507 0.531
## # ℹ 1 more variable: ACF1 <dbl>The RMSE values from the output show which model produces smaller one-step forecast errors. The model with the lower RMSE is preferred because it provides more accurate forecasts. Typically, the damped trend model is preferred if it has the smaller RMSE, since damping prevents the trend from increasing unrealistically over time while still capturing the multiplicative seasonal pattern.
d. Check that the residuals from the best method look like white
noise.best_fit <- myseries |>
model(ETS(Turnover ~ error("M") + trend("Ad") + season("M")))
best_fit |>
gg_tsresiduals()## Warning: `gg_tsresiduals()` was deprecated in feasts 0.4.2.
## ℹ Please use `ggtime::gg_tsresiduals()` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
The residual diagnostics suggest that the residuals behave like white noise. The time plot of the residuals shows values randomly scattered around zero with no clear trend or pattern. The ACF plot shows that most autocorrelations are within the significance bounds, indicating little remaining autocorrelation. The histogram is roughly symmetric and centered near zero, suggesting the residuals are approximately normally distributed. Together, these results indicate that the model has captured the main structure of the series and the remaining residuals are mostly random.
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](https://otexts.com/fpp3/toolbox-exercises.html#toolbox-exercises)?train <- myseries |>
filter(year(Month) <= 2010)
test <- myseries |>
filter(year(Month) > 2010)
fit <- train |>
model(
hw = ETS(Turnover ~ error("M") + trend("Ad") + season("M")),
snaive = SNAIVE(Turnover)
)
fc <- fit |>
forecast(new_data = test)
accuracy(fc, myseries)## # A tibble: 2 × 12
## .model State Industry .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 hw Norther… Clothin… Test 0.163 1.15 0.878 0.194 6.45 0.960 0.949 0.501
## 2 snaive Norther… Clothin… Test 0.836 1.55 1.24 5.94 9.06 1.36 1.28 0.601The test set RMSE values from accuracy show which model performs better. If the Holt Winters multiplicative model with a damped trend has a lower RMSE than the seasonal naive model then it provides more accurate forecasts. Otherwise the seasonal naive model performs better for this series.
- 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 |>
features(Turnover, guerrero) |>
pull(lambda_guerrero)
fit_stl <- train |>
model(
stl_ets = decomposition_model(
STL(box_cox(Turnover, lambda)),
ETS(season_adjust)
)
)
fc_stl <- fit_stl |>
forecast(new_data = test)
accuracy(fc_stl, myseries)## # A tibble: 1 × 12
## .model State Industry .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 stl_ets Northe… Clothin… Test -3.93 4.44 3.93 -29.8 29.8 4.30 3.66 0.879This approach stabilizes the variance of the retail turnover series using a Box Cox transformation. STL decomposition then separates the seasonal component from the data. An ETS model is fitted to the seasonally adjusted retail series and forecasts are generated. The test RMSE can then be compared with the Holt Winters and seasonal naive models to see which method produces the most accurate forecasts.