There are binary versions available but the source versions are later:
binary source needs_compilation
hardhat 1.3.1 1.4.1 FALSE
recipes 1.0.10 1.1.1 FALSE
caret 6.0-94 7.0-1 TRUE
Binaries will be installed
These will not be installed
package 'caret' successfully unpacked and MD5 sums checked
The downloaded binary packages are in
C:\Users\nbuell\AppData\Local\Temp\Rtmp8KCDAI\downloaded_packages
set.seed(9219)
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.
# Filter data for pigs in Victoriapigs <- aus_livestock |>filter(Animal =="Pigs", State =="Victoria") |>select(Count)# Fit ETS model and extract parametersfit <- pigs |>model(ETS(Count ~error("A") +trend("N") +season("N")))tidy_fit <-tidy(fit)alpha_val <- tidy_fit |>filter(term =="alpha") |>pull(estimate)l0_val <- tidy_fit |>filter(term =="l[0]") |>pull(estimate)# Generate forecasts for next 4 monthsfc <- fit |>forecast(h =4)# Plot the forecastsfc |>autoplot(pigs) +labs(y ="Number of pigs slaughtered",title ="Forecast for pig slaughter in Victoria" ) +guides(colour ="none") +theme_classic()
The optimal values of α and ℓ0, are 0.3221247 and 1.006466^{5}, respectively. I also plot the data and forecasts for the 4-month period after the data ends in December, 2018.
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.
# Get residuals standard deviations <-augment(fit) |>pull(.resid) |>sd()# Get first forecast valuefirst_forecast <- fc |>head(1) |>pull(.mean)# Calculate manual prediction intervallower <- first_forecast -1.96* supper <- first_forecast +1.96* sinterval <-c(lower, upper)interval
[1] 76871.01 113502.10
# Compare with R's intervalfc |>head(1) |>hilo() |>pull(`95%`)
<hilo[1]>
[1] [76854.79, 113518.3]95
The R interval is wider than the interval I manually calculated.
8.5
Data set global_economy contains the annual Exports from many countries. Select one country to analyse.
Plot the Exports series and discuss the main features of the data.
japanese_exports <- global_economy |>filter(Country =="Japan") |>select(Exports)japanese_exports |>ggplot(aes(Year, Exports)) +geom_line() +labs(title ="Exports from Japan generally rise, but were down significantly in the 90's.",x ="Year",y ="Exports of goods and services (% of GDP)",subtitle ="No seasonality is evident." ) +theme_classic()
Japanese exports trend up, but were down significantly in the 1990’s.
Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
First, I handle missing values.
japanese_exports |>summary()
Exports Year
Min. : 8.972 Min. :1960
1st Qu.: 9.952 1st Qu.:1974
Median :10.723 Median :1988
Mean :11.919 Mean :1988
3rd Qu.:13.893 3rd Qu.:2003
Max. :17.589 Max. :2017
NA's :1
japanese_exports |>filter(is.na(Exports) )
# A tsibble: 1 x 2 [1Y]
Exports Year
<dbl> <dbl>
1 NA 2017
clean_exports <- japanese_exports |>drop_na()
Since only the last value is missing (in 2017), I drop this missing value from the data, and forecast the next 10 years of exports.
fit <- clean_exports |>model(ETS(Exports ~error("A") +trend("N") +season("N")))fc <- fit |>forecast(h =10)fc |>autoplot(japanese_exports) +labs(title ="Forecasts of Japanese exports",x ="Year",y ="Exports of goods and services (% of GDP)" ) +theme_classic()
Compute the RMSE values for the training data.
rmse <- fit |>accuracy() |>select(RMSE)
The RMSE is 1.2633421.
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.
# A tibble: 2 × 2
.model RMSE
<chr> <dbl>
1 SES 1.26
2 Holt 1.26
The RMSE for SES (ETS(A,N,N)) is higher than for Holt’s method (ETS(A,A,N)). Holt’s method allows for trending data by adding a trend parameter, and the simpler SES model performs worse in this case. This suggests that explicitly modeling the trend improves forecast accuracy for this data.
Compare the forecasts from both methods. Which do you think is best?
# Fit both models and generate forecastsfit_both <- clean_exports |>model(SES =ETS(Exports ~error("A") +trend("N") +season("N")),Holt =ETS(Exports ~error("A") +trend("A") +season("N")) )fc_both <- fit_both |>forecast(h =10)# Generate and plot forecastsfc_both |>autoplot(japanese_exports, level =c(95), alpha =0.7) +labs(title ="Comparing SES and Holt's forecasts for Japanese exports",x ="Year",y ="Exports of goods and services (% of GDP)" ) +theme_classic()
The Holt method forecast trends up, while SES remains constant over the 10-year forecast. I think the Holt method is best because it takes the upward trend into account.
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.
# A tibble: 2 × 3
.model lower upper
<chr> <dbl> <dbl>
1 SES 13.7 18.7
2 Holt 13.9 18.8
# Get R's intervals for comparisonfc_both |>filter(Year ==min(Year)) |>hilo() |>select(.model, `95%`)
# A tsibble: 2 x 3 [1Y]
# Key: .model [2]
.model `95%` Year
<chr> <hilo> <dbl>
1 SES [13.68394, 18.72539]95 2017
2 Holt [13.80675, 18.92479]95 2017
The intervals using RMSE are wider than R’s calculated intervals.
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.]
First, I plot the training Chinese GDP data below.
# Get datachinese_gdp <- global_economy |>filter(Country =="China") |>mutate(GDP = GDP /1e12) |>select(GDP)# Plot training datachinese_gdp |>autoplot() +labs(title ="China's GDP",subtitle ="Training data",x ="Year",y ="Gross domestic product (in trillions $USD February 2019)" ) +theme_classic()
I employ the ETS statistical framework to forecast China’s GDP and investigate the following models:
I let the ETS() function select the model by minimizing the AICc.
I try damped Holt’s method.
I try using the ETS() function on the Box-Cox-transformed data.
I report on the fit and accuracy of the 20-year forecast, and plot all models below.
# Forecast 20 yearsfc <- fit |>forecast(h =10)# Get accuracyfit |>accuracy()
# A tibble: 3 × 10
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 auto Training 0.0358 0.196 0.0951 2.13 7.25 0.439 0.467 0.209
2 damped Training 0.0433 0.195 0.0933 2.84 7.53 0.430 0.465 0.0228
3 boxcox_auto Training -0.0312 0.292 0.127 0.738 7.19 0.584 0.696 0.669
# Generate and plot forecastsfc |>autoplot(chinese_gdp, level =c(95), alpha =0.7) +labs(title ="Comparing forecasts of China's GDP",x ="Year",y ="Gross domestic product (in $USD February 2019)" ) +theme_classic()
The model selected for the un-transformed data is ETS(M,A,N) and the model selected for the Box-Cox-transformed data is ETS(A,A,N). From the report of the fit, the model using the transformed data appears best (lowest absolute AIC, AICc, and BIC). It also has the best accuracy (lowest absolute RMSE and other accuracy measures). However, the damped method appears to have the best forecasts according to the RMSE of the forecasts.
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?
gas <- aus_production |>select(Gas)gas |>autoplot() +theme_classic() +labs(subtitle ="Magnitude of seasonal variation start going up in the 70's.",x ="",y ="Gas production (petajoules)" ) +theme(axis.text.x =element_text(angle =45, hjust =1)) +scale_x_yearquarter(date_breaks ="2 year", date_labels ="%Y")
gas |>model(classical_decomposition(Gas, type ="multiplicative")) |>components() |>autoplot() +labs(title ="Multiplicative decomposition shows random residuals.",x ="",y ="Gas production (petajoules)" ) +theme_classic()
gas |>model(classical_decomposition(Gas, type ="additive")) |>components() |>autoplot() +labs(title ="Additive decomposition residuals have a noticeable pattern.",subtitle ="This indicates that additive is not the right fit--it's missing something.",x ="",y ="Gas production (petajoules)" ) +theme_classic()
Since the magnitude of seasonality increases with time, there is multiplicative seasonality. The pattern of the residuals component of the additive decomposition indicate that it’s missing a pattern in the data; the residuals component of the multiplicative decomposition looks better and more random.
Next, I find an ETS model and experiment with making the trend damped.
# Fit datafit <- gas |>model(auto =ETS(Gas),damped =ETS(Gas ~error("M") +trend("Ad") +season("M")) )# Check fitreport(fit)
# Forecast 5 yearsfc <- fit |>forecast(h =4*5)# Generate and plot forecastsfc |>autoplot(gas, level =c(95), alpha =0.7) +labs(title ="Comparing forecasts of Australian gas production",x ="",y ="Gas production (petajoules)" ) +theme_classic()
The damped trend has a higher AICc than the automatically chosen ETS() model, ETS(M,A,M), indicating that it’s not as good of a fit. The plot above shows both forecasts, which appear similar. I zoom in on the most recent years below:
fc |>autoplot( gas |>filter(Quarter >yearquarter("2000 Q1")),level =c(95),alpha =0.5 ) +labs(title ="Comparing forecasts of Australian gas production from 2000 to 2015",x ="",y ="Gas production (petajoules)" ) +theme_classic()
Even after zooming in, the forecasts are very similar. You can tell that the un-damped trend has a wider CI than the damped trend further out in the forecast, and narrower CI in its first forecasts. Damping does not appear to improve the forecast.
8.8
Recall your retail time series data (from Exercise 7 in Section 2.10).
Why is multiplicative seasonality necessary for this series?
The RMSE for the damped model is higher (i.e., worse) than that of the un-damped model.
Check that the residuals from the best method look like white noise.
fit |>select(hw_mult_damped) |>gg_tsresiduals() +labs(title ="Residuals look good.",subtitle ="They appear normally distributed, random and centered around 0, and not autocorrelated." )
The innovation residuals appear randomly distributed, centered around 0. The distribution of residuals look normal. The ACF includes a few outliers, but looks OK.
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?
# Create a training dataset consisting of observations before 2011myseries_train <- turnover |>filter(year(Month) <2011)# Fit a seasonal naïve model using SNAIVE() and damped holt-winters applied to training data (myseries_train)fit <- myseries_train |>model(SNAIVE(Turnover),hw_mult_damped =ETS(Turnover ~error("M") +trend("Ad") +season("M")) )# get accuracy of model fitfit |>accuracy() |>select(.model, RMSE)
Yes, the RMSEs from the Holt-Winters’ damped model and forecasts beats (RMSE is smaller) the RSME from the seasonal naive model using the training data to the end of 2010.
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?
# Determine appropriate lambda for box cox transformationlambda <- myseries_train |>features(Turnover, features = guerrero) |>pull(lambda_guerrero)# Plot the decomposition of transformed datamyseries_train |>mutate(transformed =box_cox(Turnover, lambda)) |>model(STL(transformed ~season(window ="periodic"), robust =TRUE)) |>components() |>autoplot() +labs(title ="STL decomposition of transformed retail turnover training data") +theme_classic()
# Fit ETS model to seasonally adjusted datafit <- myseries_train |>model(# ETS of Season adj STL of box-coxseason_adj =decomposition_model(STL(box_cox(Turnover, lambda) ~season(window ="periodic")),ETS(season_adjust) ),# Holt-Winters' dampedhw_mult_damped =ETS(Turnover ~error("M") +trend("Ad") +season("M")) )# Get accuracy of model fitfit |>accuracy()
# A tibble: 2 × 10
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 season_adj Training -0.0497 3.83 2.82 -0.229 4.51 0.576 0.590 0.00426
2 hw_mult_damped Training 0.416 3.88 2.76 0.356 4.48 0.564 0.598 0.130
# Generate forecasts and get accuracy of forecastsfit |>forecast(new_data =anti_join(turnover, myseries_train)) |>accuracy(turnover)
# A tibble: 2 × 10
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 hw_mult_damped Test 13.4 19.8 16.2 9.59 12.5 3.30 3.05 0.707
2 season_adj Test -45.9 79.2 47.1 -34.0 34.9 9.62 12.2 0.325
The RMSE for the forecast of the ETS of the seasonally adjusted data from the STL decomposition of the transformed series is much higher (worse) than that of the Holt-Winters’ damped multiplicative method that was my best previous forecast of the test set (the RMSE of the model fits are very similar).
